source: git/ntl/src/xdouble.c @ 2cfffe

spielwiese
Last change on this file since 2cfffe was 2cfffe, checked in by Hans Schönemann <hannes@…>, 21 years ago
This commit was generated by cvs2svn to compensate for changes in r6316, which included commits to RCS files with non-trunk default branches. git-svn-id: file:///usr/local/Singular/svn/trunk@6317 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 14.9 KB
Line 
1
2#include <NTL/xdouble.h>
3#include <NTL/RR.h>
4
5
6#include <NTL/new.h>
7
8NTL_START_IMPL
9
10
11
12long xdouble::oprec = 10;
13
14void xdouble::SetOutputPrecision(long p)
15{
16   if (p < 1) p = 1;
17
18   if (p >= (1L << (NTL_BITS_PER_LONG-4))) 
19      Error("xdouble: output precision too big");
20
21   oprec = p;
22}
23
24void xdouble::normalize() 
25{
26   if (x == 0) 
27      e = 0;
28   else if (x > 0) {
29      while (x < NTL_XD_HBOUND_INV) { x *= NTL_XD_BOUND; e--; }
30      while (x > NTL_XD_HBOUND) { x *= NTL_XD_BOUND_INV; e++; }
31   }
32   else {
33      while (x > -NTL_XD_HBOUND_INV) { x *= NTL_XD_BOUND; e--; }
34      while (x < -NTL_XD_HBOUND) { x *= NTL_XD_BOUND_INV; e++; }
35   }
36
37   if (e >= (1L << (NTL_BITS_PER_LONG-4)))
38      Error("xdouble: overflow");
39
40   if (e <= -(1L << (NTL_BITS_PER_LONG-4)))
41      Error("xdouble: underflow");
42}
43   
44
45
46xdouble to_xdouble(double a)
47{
48   if (a == 0 || a == 1 || (a > 0 && a >= NTL_XD_HBOUND_INV && a <= NTL_XD_HBOUND)
49       || (a < 0 && a <= -NTL_XD_HBOUND_INV && a >= -NTL_XD_HBOUND)) {
50     
51      return xdouble(a, 0); 
52
53   }
54
55   if (!IsFinite(&a))
56      Error("double to xdouble conversion: non finite value");
57
58   xdouble z = xdouble(a, 0);
59   z.normalize();
60   return z;
61}
62
63
64void conv(double& xx, const xdouble& a)
65{
66   double x;
67   long e;
68
69   x = a.x;
70   e = a.e;
71
72   while (e > 0) { x *= NTL_XD_BOUND; e--; }
73   while (e < 0) { x *= NTL_XD_BOUND_INV; e++; }
74
75   xx = x;
76}
77
78
79
80
81xdouble operator+(const xdouble& a, const xdouble& b)
82{
83   xdouble z;
84
85   if (a.x == 0) 
86      return b;
87
88   if (b.x == 0)
89     return a;
90     
91
92   if (a.e == b.e) {
93      z.x = a.x + b.x;
94      z.e = a.e;
95      z.normalize();
96      return z;
97   }
98   else if (a.e > b.e) {
99      if (a.e > b.e+1)
100         return a;
101
102      z.x = a.x + b.x*NTL_XD_BOUND_INV;
103      z.e = a.e;
104      z.normalize();
105      return z;
106   }
107   else {
108      if (b.e > a.e+1)
109         return b;
110
111      z.x = a.x*NTL_XD_BOUND_INV + b.x;
112      z.e = b.e;
113      z.normalize();
114      return z;
115   }
116}
117
118
119xdouble operator-(const xdouble& a, const xdouble& b)
120{
121   xdouble z;
122
123   if (a.x == 0)
124      return -b;
125
126   if (b.x == 0)
127      return a;
128
129   if (a.e == b.e) {
130      z.x = a.x - b.x;
131      z.e = a.e;
132      z.normalize();
133      return z;
134   }
135   else if (a.e > b.e) {
136      if (a.e > b.e+1)
137         return a;
138
139      z.x = a.x - b.x*NTL_XD_BOUND_INV;
140      z.e = a.e;
141      z.normalize();
142      return z;
143   }
144   else {
145      if (b.e > a.e+1)
146         return -b;
147
148      z.x = a.x*NTL_XD_BOUND_INV - b.x;
149      z.e = b.e;
150      z.normalize();
151      return z;
152   }
153}
154
155xdouble operator-(const xdouble& a)
156{
157   xdouble z;
158   z.x = -a.x;
159   z.e = a.e;
160   return z;
161}
162
163xdouble operator*(const xdouble& a, const xdouble& b)
164{
165   xdouble z;
166
167   z.e = a.e + b.e;
168   z.x = a.x * b.x;
169   z.normalize();
170   return z;
171}
172
173xdouble operator/(const xdouble& a, const xdouble& b)
174{
175   xdouble z;
176
177   if (b.x == 0) Error("xdouble division by 0");
178
179   z.e = a.e - b.e;
180   z.x = a.x / b.x;
181   z.normalize();
182   return z;
183}
184
185
186
187long compare(const xdouble& a, const xdouble& b)
188{
189   xdouble z = a - b;
190
191   if (z.x < 0)
192      return -1;
193   else if (z.x == 0)
194      return 0;
195   else
196      return 1;
197}
198
199long sign(const xdouble& z)
200{
201   if (z.x < 0)
202      return -1;
203   else if (z.x == 0)
204      return 0;
205   else
206      return 1;
207}
208   
209
210
211xdouble trunc(const xdouble& a)
212{
213   if (a.x >= 0)
214      return floor(a);
215   else
216      return ceil(a);
217}
218
219
220xdouble floor(const xdouble& aa)
221{
222   xdouble z;
223
224   xdouble a = aa;
225   ForceToMem(&a.x);
226
227   if (a.e == 0) {
228      z.x = floor(a.x);
229      z.e = 0;
230      z.normalize();
231      return z;
232   }
233   else if (a.e > 0) {
234      return a;
235   }
236   else {
237      if (a.x < 0)
238         return to_xdouble(-1);
239      else
240         return to_xdouble(0);
241   }
242}
243
244xdouble ceil(const xdouble& aa)
245{
246   xdouble z;
247
248   xdouble a = aa;
249   ForceToMem(&a.x);
250
251   if (a.e == 0) {
252      z.x = ceil(a.x);
253      z.e = 0;
254      z.normalize();
255      return z;
256   }
257   else if (a.e > 0) {
258      return a;
259   }
260   else {
261      if (a.x < 0)
262         return to_xdouble(0);
263      else
264         return to_xdouble(1);
265   }
266}
267
268xdouble to_xdouble(const ZZ& a)
269{
270   long old_p = RR::precision();
271   RR::SetPrecision(NTL_DOUBLE_PRECISION);
272   
273   static RR t;
274   conv(t, a);
275
276   double x;
277   conv(x, t.mantissa());
278
279   xdouble y, z, res;
280
281   conv(y, x);
282   power2(z, t.exponent());
283
284   res = y*z;
285
286   RR::SetPrecision(old_p);
287
288   return res;
289}
290
291void conv(ZZ& x, const xdouble& a)
292{
293   xdouble b = floor(a);
294   long old_p = RR::precision();
295   RR::SetPrecision(NTL_DOUBLE_PRECISION);
296   static RR t;
297   conv(t, b);
298   conv(x, t);
299   RR::SetPrecision(old_p);
300}
301
302
303xdouble fabs(const xdouble& a)
304{
305   xdouble z;
306
307   z.e = a.e;
308   z.x = fabs(a.x);
309   return z;
310}
311
312xdouble sqrt(const xdouble& a)
313{
314   if (a == 0)
315      return to_xdouble(0);
316
317   if (a < 0)
318      Error("xdouble: sqrt of negative number");
319
320   xdouble t;
321
322   if (a.e & 1) {
323      t.e = (a.e - 1)/2;
324      t.x = sqrt(a.x * NTL_XD_BOUND);
325   }
326   else {
327      t.e = a.e/2;
328      t.x = sqrt(a.x);
329   }
330
331   t.normalize();
332
333   return t;
334}
335     
336
337void power(xdouble& z, const xdouble& a, const ZZ& e)
338{
339   xdouble b, res;
340
341   b = a;
342
343   res = 1;
344   long n = NumBits(e);
345   long i;
346
347   for (i = n-1; i >= 0; i--) {
348      res = res * res;
349      if (bit(e, i))
350         res = res * b;
351   }
352
353   if (sign(e) < 0) 
354      z = 1/res;
355   else
356      z = res;
357}
358
359
360
361
362void power(xdouble& z, const xdouble& a, long e)
363{
364   static ZZ E;
365   E = e;
366   power(z, a, E);
367}
368   
369
370   
371
372
373void power2(xdouble& z, long e)
374{
375   long hb = NTL_XD_HBOUND_LOG;
376   long b = 2*hb;
377
378   long q, r;
379
380   q = e/b;
381   r = e%b;
382
383   while (r >= hb) {
384      r -= b;
385      q++;
386   }
387
388   while (r < -hb) {
389      r += b;
390      q--;
391   }
392
393   if (q >= (1L << (NTL_BITS_PER_LONG-4)))
394      Error("xdouble: overflow");
395
396   if (q <= -(1L << (NTL_BITS_PER_LONG-4)))
397      Error("xdouble: underflow");
398
399   int rr = r;
400   double x = ldexp(1.0, rr);
401
402   z.x = x;
403   z.e = q;
404}
405
406
407void MulAdd(xdouble& z, const xdouble& a, const xdouble& b, const xdouble& c)
408// z = a + b*c
409{
410   double x;
411   long e;
412
413   e = b.e + c.e;
414   x = b.x * c.x;
415
416   if (x == 0) { 
417      z = a;
418      return;
419   }
420
421   if (a.x == 0) {
422      z.e = e;
423      z.x = x;
424      z.normalize();
425      return;
426   }
427     
428
429   if (a.e == e) {
430      z.x = a.x + x;
431      z.e = e;
432      z.normalize();
433      return;
434   }
435   else if (a.e > e) {
436      if (a.e > e+1) {
437         z = a;
438         return;
439      }
440
441      z.x = a.x + x*NTL_XD_BOUND_INV;
442      z.e = a.e;
443      z.normalize();
444      return;
445   }
446   else {
447      if (e > a.e+1) {
448         z.x = x;
449         z.e = e;
450         z.normalize();
451         return;
452      }
453
454      z.x = a.x*NTL_XD_BOUND_INV + x;
455      z.e = e;
456      z.normalize();
457      return;
458   }
459}
460
461void MulSub(xdouble& z, const xdouble& a, const xdouble& b, const xdouble& c)
462// z = a - b*c
463{
464   double x;
465   long e;
466
467   e = b.e + c.e;
468   x = b.x * c.x;
469
470   if (x == 0) { 
471      z = a;
472      return;
473   }
474
475   if (a.x == 0) {
476      z.e = e;
477      z.x = -x;
478      z.normalize();
479      return;
480   }
481     
482
483   if (a.e == e) {
484      z.x = a.x - x;
485      z.e = e;
486      z.normalize();
487      return;
488   }
489   else if (a.e > e) {
490      if (a.e > e+1) {
491         z = a;
492         return;
493      }
494
495      z.x = a.x - x*NTL_XD_BOUND_INV;
496      z.e = a.e;
497      z.normalize();
498      return;
499   }
500   else {
501      if (e > a.e+1) {
502         z.x = -x;
503         z.e = e;
504         z.normalize();
505         return;
506      }
507
508      z.x = a.x*NTL_XD_BOUND_INV - x;
509      z.e = e;
510      z.normalize();
511      return;
512   }
513}
514
515double log(const xdouble& a)
516{
517   static double LogBound = log(NTL_XD_BOUND);
518   if (a.x <= 0) {
519      Error("log(xdouble): argument must be positive");
520   }
521
522   return log(a.x) + a.e*LogBound;
523}
524
525xdouble xexp(double x)
526{
527   const double LogBound = log(NTL_XD_BOUND);
528
529   double y = x/LogBound;
530   double iy = floor(y+0.5);
531
532   if (iy >= (1L << (NTL_BITS_PER_LONG-4)))
533      Error("xdouble: overflow");
534
535   if (iy <= -(1L << (NTL_BITS_PER_LONG-4)))
536      Error("xdouble: underflow");
537
538
539   double fy = y - iy;
540
541   xdouble res;
542   res.e = long(iy);
543   res.x = exp(fy*LogBound);
544   res.normalize();
545   return res;
546}
547
548/**************  input / output routines **************/
549
550
551void ComputeLn2(RR&);
552void ComputeLn10(RR&);
553
554long ComputeMax10Power()
555{
556   long old_p = RR::precision();
557   RR::SetPrecision(NTL_BITS_PER_LONG);
558
559   RR ln2, ln10;
560   ComputeLn2(ln2);
561   ComputeLn10(ln10);
562
563   long k = to_long( to_RR(1L << (NTL_BITS_PER_LONG-5)) * ln2 / ln10 );
564
565   RR::SetPrecision(old_p);
566   return k;
567}
568
569
570xdouble PowerOf10(const ZZ& e)
571{
572   static long init = 0;
573   static xdouble v10k;
574   static long k;
575
576   if (!init) {
577      long old_p = RR::precision();
578      k = ComputeMax10Power();
579      RR::SetPrecision(NTL_DOUBLE_PRECISION);
580      v10k = to_xdouble(power(to_RR(10), k)); 
581      RR::SetPrecision(old_p);
582      init = 1;
583   }
584
585   ZZ e1;
586   long neg;
587
588   if (e < 0) {
589      e1 = -e;
590      neg = 1;
591   }
592   else {
593      e1 = e;
594      neg = 0;
595   }
596
597   long r;
598   ZZ q;
599
600   r = DivRem(q, e1, k);
601
602   long old_p = RR::precision();
603   RR::SetPrecision(NTL_DOUBLE_PRECISION);
604   xdouble x1 = to_xdouble(power(to_RR(10), r));
605   RR::SetPrecision(old_p);
606
607   xdouble x2 = power(v10k, q);
608   xdouble x3 = x1*x2;
609
610   if (neg) x3 = 1/x3;
611
612   return x3;
613}
614
615
616
617
618ostream& operator<<(ostream& s, const xdouble& a)
619{
620   if (a == 0) {
621      s << "0";
622      return s;
623   }
624
625   long old_p = RR::precision();
626   long temp_p = long(log(fabs(log(fabs(a))) + 1.0)/log(2.0)) + 10; 
627
628   RR::SetPrecision(temp_p);
629
630   RR ln2, ln10, log_2_10;
631   ComputeLn2(ln2);
632   ComputeLn10(ln10);
633   log_2_10 = ln10/ln2;
634   ZZ log_10_a = to_ZZ(
635  (to_RR(a.e)*to_RR(2*NTL_XD_HBOUND_LOG) + log(fabs(a.x))/log(2.0))/log_2_10);
636
637   RR::SetPrecision(old_p);
638
639   xdouble b;
640   long neg;
641
642   if (a < 0) {
643      b = -a;
644      neg = 1;
645   }
646   else {
647      b = a;
648      neg = 0;
649   }
650
651   ZZ k = xdouble::OutputPrecision() - log_10_a;
652
653   xdouble c, d;
654
655   c = PowerOf10(to_ZZ(xdouble::OutputPrecision()));
656   d = PowerOf10(log_10_a);
657
658   b = b / d;
659   b = b * c;
660
661   while (b < c) {
662      b = b * 10.0;
663      k++;
664   }
665
666   while (b >= c) {
667      b = b / 10.0;
668      k--;
669   }
670
671   b = b + 0.5;
672   k = -k;
673
674   ZZ B;
675   conv(B, b);
676
677   long bp_len = xdouble::OutputPrecision()+10;
678
679   char *bp = NTL_NEW_OP char[bp_len];
680
681   if (!bp) Error("xdouble output: out of memory");
682
683   long len, i;
684
685   len = 0;
686   do {
687      if (len >= bp_len) Error("xdouble output: buffer overflow");
688      bp[len] = DivRem(B, B, 10) + '0';
689      len++;
690   } while (B > 0);
691
692   for (i = 0; i < len/2; i++) {
693      char tmp;
694      tmp = bp[i];
695      bp[i] = bp[len-1-i];
696      bp[len-1-i] = tmp;
697   }
698
699   i = len-1;
700   while (bp[i] == '0') i--;
701
702   k += (len-1-i);
703   len = i+1;
704
705   bp[len] = '\0';
706
707   if (k > 3 || k < -len - 3) {
708      // use scientific notation
709
710      if (neg) s << "-";
711      s << "0." << bp << "e" << (k + len);
712   }
713   else {
714      long kk = to_long(k);
715
716      if (kk >= 0) {
717         if (neg) s << "-";
718         s << bp;
719         for (i = 0; i < kk; i++) 
720            s << "0";
721      }
722      else if (kk <= -len) {
723         if (neg) s << "-";
724         s << "0.";
725         for (i = 0; i < -len-kk; i++)
726            s << "0";
727         s << bp;
728      }
729      else {
730         if (neg) s << "-";
731         for (i = 0; i < len+kk; i++)
732            s << bp[i];
733   
734         s << ".";
735   
736         for (i = len+kk; i < len; i++)
737            s << bp[i];
738      }
739   }
740
741   delete [] bp;
742   return s;
743}
744
745istream& operator>>(istream& s, xdouble& x)
746{
747   long c;
748   long sign;
749   ZZ a, b;
750
751   if (!s) Error("bad xdouble input");
752
753   c = s.peek();
754   while (c == ' ' || c == '\n' || c == '\t') {
755      s.get();
756      c = s.peek();
757   }
758
759   if (c == '-') {
760      sign = -1;
761      s.get();
762      c = s.peek();
763   }
764   else
765      sign = 1;
766
767   long got1 = 0;
768   long got_dot = 0;
769   long got2 = 0;
770
771   a = 0;
772   b = 1;
773
774   if (c >= '0' && c <= '9') {
775      got1 = 1;
776
777      while (c >= '0' && c <= '9') {
778         mul(a, a, 10);
779         add(a, a, c-'0');
780         s.get();
781         c = s.peek();
782      }
783   }
784
785   if (c == '.') {
786      got_dot = 1;
787
788      s.get();
789      c = s.peek();
790
791      if (c >= '0' && c <= '9') {
792         got2 = 1;
793   
794         while (c >= '0' && c <= '9') {
795            mul(a, a, 10);
796            add(a, a, c-'0');
797            mul(b, b, 10);
798            s.get();
799            c = s.peek();
800         }
801      }
802   }
803
804   if (got_dot && !got1 && !got2)  Error("bad xdouble input");
805
806   ZZ e;
807
808   long got_e = 0;
809   long e_sign;
810
811   if (c == 'e' || c == 'E') {
812      got_e = 1;
813
814      s.get();
815      c = s.peek();
816
817      if (c == '-') {
818         e_sign = -1;
819         s.get();
820         c = s.peek();
821      }
822      else if (c == '+') {
823         e_sign = 1;
824         s.get();
825         c = s.peek();
826      }
827      else
828         e_sign = 1;
829
830      if (c < '0' || c > '9') Error("bad xdouble input");
831
832      e = 0;
833      while (c >= '0' && c <= '9') {
834         mul(e, e, 10);
835         add(e, e, c-'0');
836         s.get();
837         c = s.peek();
838      }
839   }
840
841   if (!got1 && !got2 && !got_e) Error("bad xdouble input");
842
843   xdouble t1, t2, v;
844
845   if (got1 || got2) {
846      conv(t1, a);
847      conv(t2, b);
848      v = t1/t2;
849   }
850   else
851      v = 1;
852
853   if (sign < 0)
854      v = -v;
855
856   if (got_e) {
857      if (e_sign < 0) negate(e, e);
858      t1 = PowerOf10(e);
859      v = v * t1;
860   }
861
862   x = v;
863   return s;
864}
865
866
867xdouble to_xdouble(const char *s)
868{
869   long c;
870   long sign;
871   ZZ a, b;
872   long i=0;
873
874   if (!s) Error("bad xdouble input");
875
876   c = s[i];
877   while (c == ' ' || c == '\n' || c == '\t') {
878      i++;
879      c = s[i];
880   }
881
882   if (c == '-') {
883      sign = -1;
884      i++;
885      c = s[i];
886   }
887   else
888      sign = 1;
889
890   long got1 = 0;
891   long got_dot = 0;
892   long got2 = 0;
893
894   a = 0;
895   b = 1;
896
897   if (c >= '0' && c <= '9') {
898      got1 = 1;
899
900      while (c >= '0' && c <= '9') {
901         mul(a, a, 10);
902         add(a, a, c-'0');
903         i++;
904         c = s[i];
905      }
906   }
907
908   if (c == '.') {
909      got_dot = 1;
910
911      i++;
912      c = s[i];
913
914      if (c >= '0' && c <= '9') {
915         got2 = 1;
916   
917         while (c >= '0' && c <= '9') {
918            mul(a, a, 10);
919            add(a, a, c-'0');
920            mul(b, b, 10);
921            i++;
922            c = s[i];
923         }
924      }
925   }
926
927   if (got_dot && !got1 && !got2)  Error("bad xdouble input");
928
929   ZZ e;
930
931   long got_e = 0;
932   long e_sign;
933
934   if (c == 'e' || c == 'E') {
935      got_e = 1;
936
937      i++;
938      c = s[i];
939
940      if (c == '-') {
941         e_sign = -1;
942         i++;
943         c = s[i];
944      }
945      else if (c == '+') {
946         e_sign = 1;
947         i++;
948         c = s[i];
949      }
950      else
951         e_sign = 1;
952
953      if (c < '0' || c > '9') Error("bad xdouble input");
954
955      e = 0;
956      while (c >= '0' && c <= '9') {
957         mul(e, e, 10);
958         add(e, e, c-'0');
959         i++;
960         c = s[i];
961      }
962   }
963
964   if (!got1 && !got2 && !got_e) Error("bad xdouble input");
965
966   xdouble t1, t2, v;
967
968   if (got1 || got2) {
969      conv(t1, a);
970      conv(t2, b);
971      v = t1/t2;
972   }
973   else
974      v = 1;
975
976   if (sign < 0)
977      v = -v;
978
979   if (got_e) {
980      if (e_sign < 0) negate(e, e);
981      t1 = PowerOf10(e);
982      v = v * t1;
983   }
984
985   return v;
986}
987
988NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.