source: trunk/hiptop/RetroStatus/dgs/libs/hipfloat/hipfloat.java@ 207

Last change on this file since 207 was 207, checked in by Nicholas Riley, 18 years ago

rpn: Brian Swetland's RPN calculator using Dan Sachs's hipfloat
library. Changes for hiptop OS 2.3; some minor UI cleanups.
File size: 15.7 KB
Line 
1package dgs.libs.hipfloat;
2
3// Copyright 2003, Daniel Grobe Sachs. All Rights Reserved.
4// See LICENSE for redistribution terms
5//
6// Some algorithms borrowed from GNU BC, but all code was rewritten.
7//
8// arctrig functions contributed by Greg Vander Rhodes <greg@vanderrhodes.com>
9
10public class hipfloat implements Comparable {
11 protected int mantissa;
12 protected int exponent;
13
14 static public final hipfloat ZERO = new hipfloat(0);
15 static public final hipfloat ONE = new hipfloat(1);
16 static public final hipfloat HALF = new hipfloat(5,-1);
17 static public final hipfloat TWO = new hipfloat(2);
18 static public final hipfloat PI = new hipfloat(314159265,-8);
19 static public final hipfloat PI2 = new hipfloat(628318530,-8);
20 static public final hipfloat HALFPI = PI.div(TWO);
21 static public final hipfloat E = new hipfloat(271828182,-8);
22 static public final hipfloat MAXEXP = new hipfloat(23026);
23
24 static public final hipfloatBadNum OVF = hipfloatBadNum.OVF;
25 static public final hipfloatBadNum NAN = hipfloatBadNum.NAN;
26 static public final hipfloatBadNum NoError = hipfloatBadNum.NoError;
27
28 static private final int fact_db[] =
29 { 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916000 };
30
31 static public final int MAX_EXP = 9999;
32
33 public hipfloat(hipfloat in)
34 {
35 hipfloat out;
36
37 mantissa = in.mantissa;
38 exponent = in.exponent;
39
40 out = normalize(this);
41
42 if( out != this ) // must be an error
43 throw new hipfloatError("Overflow",OVF);
44 }
45
46 public hipfloat(int in)
47 {
48 mantissa = in;
49 exponent = 0;
50
51 normalize(this);
52 }
53
54 public hipfloat(int man, int exp)
55 {
56 hipfloat out;
57
58 mantissa = man;
59 exponent = exp;
60
61 out = normalize(this);
62
63 if( out != this ) // must be an error
64 throw new hipfloatError("Overflow",OVF);
65 }
66
67 public hipfloat(long man, int exp)
68 {
69 while ( (exp < -(MAX_EXP+8)) && (man != 0))
70 {
71 exp += 1;
72 man /= 10;
73 }
74
75 if( man == 0 )
76 {
77 mantissa = 0;
78 exponent = -8;
79 return;
80 }
81
82 while( (Math.abs(man) >= 1000000000) && (exp < (MAX_EXP-8)) )
83 {
84 if( (Math.abs(man) < 10000000000L) && ((Math.abs(man)%10) >= 5) )
85 man += 10; // need to round up
86
87 man /= 10;
88 exp += 1;
89 }
90
91 while( (Math.abs(man) < 100000000) && (exp > -(MAX_EXP+8)) )
92 {
93 man *= 10;
94 exp -= 1;
95 }
96
97 if( (exp > (MAX_EXP-8)) || (man >= 1000000000) )
98 throw new hipfloatError("Overflow",OVF);
99
100 mantissa = (int)man;
101 exponent = exp;
102 }
103
104 public hipfloat(String in)
105 {
106 hipfloat work = fromString(in);
107
108 if( work.isError() )
109 {
110 if( work == OVF )
111 throw new hipfloatError("Overflow on conversion",
112 (hipfloatBadNum)work);
113 else if( work == NAN )
114 throw new hipfloatError("Invalid number conversion",
115 (hipfloatBadNum)work);
116 else
117 throw new hipfloatError("Unknown error converting string",
118 (hipfloatBadNum)work);
119 }
120
121 this.mantissa = work.mantissa;
122 this.exponent = work.exponent;
123 }
124
125 public int mantissa()
126 {
127 return mantissa;
128 }
129
130 public int exponent()
131 {
132 return exponent;
133 }
134
135 public hipfloat floor()
136 {
137 if( exponent < -10 ) return ZERO;
138 if( exponent > 0 ) return this;
139
140 hipfloat out = new hipfloat(this);
141
142 while( out.exponent < 0 )
143 {
144 out.exponent++;
145 out.mantissa /= 10;
146 }
147
148 return normalize(out);
149 }
150
151 public hipfloat ceil()
152 {
153 if( this.compareTo(this.floor()) == 0 )
154 return this;
155
156 return( (this.add(ONE)).floor() );
157 }
158
159 public hipfloat round()
160 {
161 return (this.add(HALF)).floor();
162 }
163
164 public int toint()
165 {
166 if( (this.exponent > 0) ||
167 ((this.exponent == 1) && (this.mantissa > 214748364)) )
168 throw new hipfloatError("Overflow",OVF);
169
170 if( this.exponent == 1 )
171 return this.mantissa * 10;
172
173 int exp = this.exponent, man = this.mantissa;
174
175 while( exp < 0 )
176 {
177 exp ++;
178 man /= 10;
179 }
180
181 return man;
182 }
183
184
185 public hipfloat add(hipfloat a)
186 {
187 if( a.isError() ) return a;
188
189 if( mantissa == 0 )
190 return a;
191
192 if( a.mantissa == 0 )
193 return this;
194
195 if( a.exponent > exponent )
196 return a.add(this);
197 else
198 {
199 int aexp, aman;
200
201 aexp = a.exponent;
202 aman = a.mantissa;
203
204 while( aexp < exponent-1 )
205 {
206 aexp++;
207 aman /= 10;
208 }
209
210 if( aexp == exponent-1 )
211 {
212 if( Math.abs(aman)%10 >= 5 )
213 aman += 10;
214
215 aman /= 10;
216 aexp++;
217 }
218
219 try
220 { return new hipfloat(mantissa+aman, exponent); }
221 catch( hipfloatError error )
222 { return error.actual_return; }
223 }
224
225 }
226
227 public hipfloat mul(hipfloat a)
228 {
229 if( a.isError() ) return a;
230
231 long man;
232 int exp;
233
234 exp = exponent + a.exponent;
235 man = (long)mantissa * a.mantissa;
236
237 try
238 { return new hipfloat(man, exp); }
239 catch( hipfloatError error )
240 { return error.actual_return; }
241 }
242
243 public hipfloat sub(hipfloat a)
244 {
245 if( a.isError() ) return a;
246
247 try
248 { return this.add(a.neg()); }
249 catch( hipfloatError error )
250 { return error.actual_return; }
251 }
252
253 public hipfloat div(hipfloat a)
254 {
255 if( a.isError() ) return a;
256
257 if( a.compareTo(ZERO) == 0 )
258 return NAN;
259
260 long divs = 1000000000000000000L/a.mantissa;;
261 int exp = -a.exponent - 18;
262
263 hipfloat b;
264
265 try
266 { b = new hipfloat(divs, exp); }
267 catch( hipfloatError error )
268 { return error.actual_return; }
269
270 return this.mul(b);
271 }
272
273 public hipfloat pow(hipfloat a)
274 {
275 if( a.isError() ) return a;
276
277 hipfloat x;
278 int i, t;
279
280 x = this;
281
282 if( a.compareTo(a.floor()) == 0 ) // if integer
283 try
284 {
285 t = (a.abs()).toint();
286
287 hipfloat out=ONE,z=x;
288
289 while( t > 0 )
290 {
291 if( 1 == (t % 2) )
292 out = out.mul(z);
293
294 z = z.mul(z);
295 t = t / 2;
296 }
297
298 return out;
299 }
300 catch( hipfloatError error )
301 {
302 // fall back to using logs
303 ;
304 }
305
306 if( x.compareTo(ZERO) < 0 )
307 return NAN; // if a isn't integer, this isn't valid
308
309 return ((x.ln()).mul(a)).exp();
310 }
311
312 public hipfloat sqrt()
313 {
314 hipfloat guess;
315 hipfloat check;
316
317 if( this.compareTo(ZERO) < 0 )
318 return NAN;
319
320 if( this.compareTo(ZERO) == 0 )
321 return ZERO;
322
323 try {
324 int i;
325
326 //if( this.compareTo(new hipfloat(1)) > 0 )
327 guess = new hipfloat(1,(exponent+8)/2);
328 //else
329 // guess = new hipfloat(1,0);
330
331 for( i = 0; i < 8; i++ )
332 {
333 check = this.div(guess);
334 guess = (check.add(guess)).mul(HALF);
335 }
336 }
337 catch( hipfloatError error )
338 { return error.actual_return; }
339
340 return guess;
341 }
342
343 // Use the Taylor series:
344 // 3 5
345 // x-1 1 x-1 1 x-1
346 // ln(x) =2 --- + - --- + - --- + ...
347 // x+1 3 x+1 5 x+1
348
349 public hipfloat ln()
350 {
351 if( this.compareTo(ZERO) < 0 )
352 return NAN;
353
354 hipfloat in = new hipfloat(this);
355 int scale = 2;
356 int i;
357
358 // bring input to the range 0.5 - 2, exclusive
359
360 while( (in.compareTo(TWO) >= 0) || (in.compareTo(HALF) <= 0) )
361 {
362 scale *= 2;
363 in = in.sqrt();
364 }
365
366 hipfloat t;
367 hipfloat t2;
368
369 hipfloat out = ZERO;
370
371 t = (in.sub(ONE)).div(in.add(ONE));
372 t2 = t.mul(t);
373
374
375 for( i = 0; i < 20; i++ )
376 {
377 out = out.add(t.div(new hipfloat(1+2*i)));
378 t = t.mul(t2);
379
380 if( t.nearTo(ZERO,10) )
381 break;
382 }
383
384 return out.mul(new hipfloat(scale));
385 }
386
387 public hipfloat exp()
388 {
389 hipfloat x = new hipfloat(this);
390 hipfloat xx = new hipfloat(0);
391 hipfloat out = new hipfloat(0);
392 int count = 0, i;
393
394 // don't waste time doing computation if it'll overflow
395
396 if( x.compareTo(MAXEXP) > 0 )
397 return OVF;
398
399 // Scale to range where we can use the Taylor series
400
401 while( x.compareTo(ONE) > 0 )
402 {
403 x = x.div(TWO);
404 count ++;
405 }
406
407 // Taylor series of 1+x+x^2/2!+x^3/3! ...
408
409 out = (new hipfloat(1)).add(x);
410 xx = x.mul(x);
411
412 for( i = 2; i < 12; i++ )
413 {
414 out = out.add(xx.div(this.factorial(i)));
415 xx = xx.mul(x);
416 }
417
418 while( count-- > 0 )
419 out = out.mul(out);
420
421 return out;
422 }
423
424 public hipfloat sin()
425 {
426 // series is x - x^3/3! + x^5/5! - ...
427
428 int i, sign = -1;
429 hipfloat x = this;
430
431 if( x.exponent > 0 )
432 return OVF;
433
434 while( x.abs().compareTo(PI) > 0 )
435 {
436 int exp = x.exponent + 7;
437 if( exp < 0 ) exp = 0;
438
439 if( x.compareTo(ZERO) > 0 )
440 x = x.sub(PI2.mul(new hipfloat(1,exp)));
441 else
442 x = x.add(PI2.mul(new hipfloat(1,exp)));
443 }
444
445 if( x.abs().compareTo(HALFPI) > 0 )
446 if( x.compareTo(ZERO) > 0 )
447 x = PI.sub(x);
448 else
449 x = PI.neg().add(x);
450
451 hipfloat out = x;
452 hipfloat f = x.mul(x);
453 hipfloat d = ONE;
454
455 for( i = 3; i < 13; i += 2)
456 {
457 x = x.mul(f);
458 d = d.mul(new hipfloat(i-1)).mul(new hipfloat(i));
459
460 if( sign == 1 )
461 out = out.add(x.div(d));
462 else
463 out = out.sub(x.div(d));
464
465 sign = -sign;
466 }
467
468 return out;
469 }
470
471 public hipfloat cos()
472 {
473 return HALFPI.sub(this).sin();
474 }
475
476 public hipfloat tan()
477 {
478 return this.sin().div(this.cos());
479 }
480
481 public hipfloat arctan ()
482 {
483 int i, sign = -1;
484 hipfloat x = this.abs ();
485 hipfloat xsign = new hipfloat (1);
486 if (this.compareTo (ZERO) < 0)
487 {
488 xsign = xsign.neg ();
489 }
490 boolean over1flag = false;
491 boolean oneflag = false;
492 if (x.compareTo (ONE) > 0)
493 {
494 x = ONE.div (x);
495 over1flag = true;
496 }
497
498 if (x.compareTo (ONE) == 0)
499 {
500 oneflag = true;
501 }
502 // use series
503 // arctan(y) = y - y^3/3 + y^5/5 - y^7/7 + y^9/9 - ...
504
505 hipfloat out = x;
506 if (!oneflag)
507 {
508
509 hipfloat f = x.mul (x);
510 hipfloat d = ONE;
511
512 for (i = 3; i < 23; i += 2)
513 {
514 x = x.mul (f);
515 d = new hipfloat (i);
516
517 if (sign == 1)
518 out = out.add (x.div (d));
519 else
520 out = out.sub (x.div (d));
521
522 sign = -sign;
523 }
524
525 if (over1flag)
526 {
527 out = PI.div (new hipfloat (2)).sub (out);
528 }
529 }
530 else
531 {
532 out = new hipfloat (785398163, -9);
533 }
534
535 return out.mul (xsign);
536 }
537
538 public hipfloat arccos ()
539 {
540 if (this.abs ().compareTo (ONE) > 0)
541 return NAN;
542
543 hipfloat x = this.abs ();
544 boolean negflag = false;
545 boolean zeroflag = false;
546 boolean oneflag = false;
547 if (this.compareTo (ZERO) < 0)
548 negflag = true;
549
550 if (this.compareTo (ZERO) == 0)
551 zeroflag = true;
552
553 if (x.compareTo (ONE) == 0)
554 oneflag = true;
555
556 hipfloat out = new hipfloat (1);
557
558 if (zeroflag)
559 out = PI.div (TWO);
560 else
561 {
562 if (oneflag)
563 out = ZERO;
564 else
565 {
566 out = out.sub (x.mul (x)).sqrt ().div (x);
567 out = out.arctan ();
568 }
569 }
570
571 if (negflag)
572 out = PI.sub (out);
573
574 return out;
575 }
576
577 public hipfloat arcsin ()
578 {
579 if (this.abs ().compareTo (ONE) > 0)
580 return NAN;
581
582 hipfloat x = this.abs ();
583 boolean negflag = false;
584 boolean zeroflag = false;
585 boolean oneflag = false;
586 if (this.compareTo (ZERO) < 0)
587 negflag = true;
588
589 if (this.compareTo (ZERO) == 0)
590 zeroflag = true;
591
592 if (x.compareTo (ONE) == 0)
593 oneflag = true;
594
595 hipfloat out = new hipfloat (1);
596
597 if (zeroflag)
598 out = ZERO;
599 else
600 {
601 if (oneflag)
602 out = PI.div (TWO);
603 else
604 {
605 out = x.div (out.sub (x.mul (x)).sqrt ());
606 out = out.arctan ();
607 }
608 }
609
610 if (negflag)
611 out = out.neg ();
612
613 return out;
614 }
615
616 public hipfloat neg()
617 {
618 return new hipfloat(-mantissa,exponent);
619 }
620
621 static public hipfloat factorial(int i)
622 {
623 if( i < 12 )
624 return new hipfloat(fact_db[i]);
625 else
626 return (new hipfloat(i)).factorial();
627 }
628
629 public hipfloat factorial()
630 {
631 if( this.compareTo(ZERO) <= 0 )
632 return new hipfloat( 1 );
633
634 return this.mul((this.sub(ONE)).factorial());
635 }
636
637 public hipfloat abs()
638 {
639 return new hipfloat(Math.abs(mantissa),exponent);
640 }
641
642 public boolean nearTo(hipfloat a, int ulp)
643 {
644 if( a.isError() )
645 return false;
646
647 int x;
648
649 hipfloat c1 = new hipfloat(this);
650 hipfloat c2 = new hipfloat(a);
651
652 while( c1.exponent < c2.exponent )
653 {
654 c2.exponent --;
655 c2.mantissa /= 10;
656 }
657
658 while( c2.exponent < c1.exponent )
659 {
660 c1.exponent --;
661 c1.mantissa /= 10;
662 }
663
664 x = c1.mantissa - c2.mantissa;
665
666 if( Math.abs(x) <= ulp )
667 return true;
668 else
669 return false;
670 }
671
672 public boolean isError()
673 {
674 return false;
675 }
676
677 public int compareTo(Object o)
678 {
679 if( ((hipfloat)(o)).isError() )
680 return 1;
681
682 hipfloat a = (hipfloat)o;
683
684 if( a.mantissa == 0 )
685 return this.mantissa;
686
687 return (this.sub(a)).mantissa;
688 }
689
690 static int returnError(char[] out, String error)
691 {
692 int i;
693
694 for( i = 0; i < error.length(); i++ )
695 out[i] = error.charAt(i);
696
697 return error.length();
698 }
699
700 public int toCharArray(char[] out, boolean Scientific)
701 {
702 int unit;
703 int frac;
704 int vexp = exponent + 8;
705 int i, j, sigfig;
706
707 char c;
708 int len = 0;
709
710 if( this == NAN )
711 return returnError(out, "Not a number");
712
713 if( this == OVF )
714 return returnError(out, "Overflow");
715
716 if( mantissa < 0 )
717 out[len++] = '-';
718
719 frac = Math.abs(mantissa);
720
721 if( !Scientific )
722 { // Can we print this number?
723 i = 100000000;
724 sigfig = 1;
725
726 while( (i > 0) && (0!=(mantissa % i)) )
727 {
728 sigfig ++;
729 i /= 10;
730 }
731
732 if( (vexp == -1) && (sigfig == 9) )
733 {
734 out[len++] = '0';
735 out[len++] = '.';
736 vexp = 8;
737 }
738
739 if( (vexp <= 8) && (vexp >= -9+sigfig))
740 { // best not scientific
741 i = 100000000;
742
743 while( vexp < 0 )
744 {
745 frac = frac/10;
746 vexp++;
747 sigfig++;
748 }
749
750 for( j = 0; j < sigfig; j++ )
751 {
752 if( vexp == -1 )
753 out[len++] = '.';
754
755 unit = frac/i;
756 frac = frac%i;
757 out[len++] = (char) ('0' + (char)unit);
758
759 vexp--;
760 i /= 10;
761 }
762
763 while( vexp >= 0 )
764 {
765 out[len++] = '0';
766 vexp--;
767 }
768
769 return len;
770 }
771 }
772
773 for( i = 100000000; i > 0 ; i /= 10 )
774 {
775 unit = frac/i;
776 frac = frac%i;
777
778 out[len++] = (char) ('0'+(char)unit);
779 if( i == 100000000 )
780 out[len++] = '.';
781 }
782
783 out[len++] = 'e';
784
785 if( vexp < 0 )
786 {
787 out[len++] = '-';
788 vexp = -vexp;
789 }
790 else
791 out[len++] = '+';
792
793
794 boolean printed = false;
795
796 for( i = 10000; i > 0; i /= 10 )
797 {
798 if( (vexp >= i) || (i == 1) || printed )
799 {
800 printed = true;
801 out[len++] = (char) ('0' + (char) (vexp/i));
802 }
803
804 vexp %= i;
805
806 }
807
808
809 return len;
810 }
811
812 public String toString(boolean Scientific)
813 {
814 char out[] = new char[64];
815 int len;
816
817 len = this.toCharArray(out, Scientific);
818
819 return new String(out,0,len);
820 }
821
822 public String toString()
823 {
824 return toString(false);
825 }
826
827 public static hipfloat fromString(String in)
828 {
829 int i, j, state = 0, exp = 0, dot = 0;
830 boolean neg = false, expneg = false;
831
832 hipfloat work = new hipfloat(0);
833
834 hipfloatBadNum error = NoError;
835
836 char c;
837
838 in = in.trim();
839 i = 0;
840
841 while( (error == NoError) && (i < in.length()) )
842 {
843 c = in.charAt(i++);
844
845 switch( state )
846 {
847 case 0: state = 1;
848
849 if( c == '-' )
850 {
851 neg = true;
852 break;
853 }
854
855 if( c == '+' )
856 break;
857
858 case 1: if( (c >= '0') && (c <= '9') )
859 {
860 if( dot == 0 )
861 work = (work.mul(new hipfloat(10))).
862 add(new hipfloat(c-'0'));
863 else
864 {
865 work = work.add(new hipfloat((c-'0'),dot));
866 dot --;
867 }
868 }
869 else if( (c == '.') && (dot == 0) )
870 dot = -1;
871 else if( (c == 'e') || (c == 'E') )
872 state = 2;
873 else
874 error = NAN;
875
876 break;
877
878 case 2: state = 3;
879
880 if( c == '-' )
881 {
882 expneg = true;
883 break;
884 }
885
886 if( c == '+' )
887 break;
888
889 case 3: if( (c >= '0') && (c <= '9') )
890 {
891 exp = exp * 10 + (c - '0');
892
893 if( exp > 10000 )
894 error = OVF;
895 }
896 else
897 error = NAN;
898 }
899 }
900
901 if( error != NoError )
902 return error;
903
904 if( neg )
905 work = work.neg();
906
907 if( expneg )
908 exp = -exp;
909
910 if( exp != 0 )
911 work = work.mul(new hipfloat(1,exp));
912
913 //System.err.println(work);
914
915 return work;
916 }
917
918 private static hipfloat normalize(hipfloat in)
919 {
920 while ( (in.exponent < -(MAX_EXP+8)) && (in.mantissa != 0))
921 {
922 in.exponent += 1;
923 in.mantissa /= 10;
924 }
925
926 if( in.mantissa == 0 )
927 {
928 in.exponent = -8;
929 return in;
930 }
931
932 while( (Math.abs(in.mantissa) >= 1000000000) && (in.exponent < (MAX_EXP-8)) )
933 {
934 if( (Math.abs(in.mantissa % 10)) >= 5 )
935 in.mantissa += 10; // need to round up
936
937 in.mantissa /= 10;
938 in.exponent += 1;
939 }
940
941 while( (Math.abs(in.mantissa) < 100000000) && (in.exponent > -(MAX_EXP+8)) )
942 {
943 in.mantissa *= 10;
944 in.exponent -= 1;
945 }
946
947 if( (in.exponent > (MAX_EXP-8)) || (in.mantissa >= 1000000000) )
948 return OVF;
949 else
950 return in;
951 }
952}
953
Note: See TracBrowser for help on using the repository browser.