source: git/ntl/src/xdouble.c @ 33a041

spielwiese
Last change on this file since 33a041 was 311902, checked in by Hans Schönemann <hannes@…>, 20 years ago
*hannes: 5.3.2-C++-fix git-svn-id: file:///usr/local/Singular/svn/trunk@7474 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 10.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 (NTL_OVERFLOW(p, 1, 0)) 
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 >= NTL_OVFBND)
38      Error("xdouble: overflow");
39
40   if (e <= -NTL_OVFBND)
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 >= NTL_OVFBND)
394      Error("xdouble: overflow");
395
396   if (q <= -NTL_OVFBND)
397      Error("xdouble: underflow");
398
399   double x = _ntl_ldexp(1.0, r);
400
401   z.x = x;
402   z.e = q;
403}
404
405
406void MulAdd(xdouble& z, const xdouble& a, const xdouble& b, const xdouble& c)
407// z = a + b*c
408{
409   double x;
410   long e;
411
412   e = b.e + c.e;
413   x = b.x * c.x;
414
415   if (x == 0) { 
416      z = a;
417      return;
418   }
419
420   if (a.x == 0) {
421      z.e = e;
422      z.x = x;
423      z.normalize();
424      return;
425   }
426     
427
428   if (a.e == e) {
429      z.x = a.x + x;
430      z.e = e;
431      z.normalize();
432      return;
433   }
434   else if (a.e > e) {
435      if (a.e > e+1) {
436         z = a;
437         return;
438      }
439
440      z.x = a.x + x*NTL_XD_BOUND_INV;
441      z.e = a.e;
442      z.normalize();
443      return;
444   }
445   else {
446      if (e > a.e+1) {
447         z.x = x;
448         z.e = e;
449         z.normalize();
450         return;
451      }
452
453      z.x = a.x*NTL_XD_BOUND_INV + x;
454      z.e = e;
455      z.normalize();
456      return;
457   }
458}
459
460void MulSub(xdouble& z, const xdouble& a, const xdouble& b, const xdouble& c)
461// z = a - b*c
462{
463   double x;
464   long e;
465
466   e = b.e + c.e;
467   x = b.x * c.x;
468
469   if (x == 0) { 
470      z = a;
471      return;
472   }
473
474   if (a.x == 0) {
475      z.e = e;
476      z.x = -x;
477      z.normalize();
478      return;
479   }
480     
481
482   if (a.e == e) {
483      z.x = a.x - x;
484      z.e = e;
485      z.normalize();
486      return;
487   }
488   else if (a.e > e) {
489      if (a.e > e+1) {
490         z = a;
491         return;
492      }
493
494      z.x = a.x - x*NTL_XD_BOUND_INV;
495      z.e = a.e;
496      z.normalize();
497      return;
498   }
499   else {
500      if (e > a.e+1) {
501         z.x = -x;
502         z.e = e;
503         z.normalize();
504         return;
505      }
506
507      z.x = a.x*NTL_XD_BOUND_INV - x;
508      z.e = e;
509      z.normalize();
510      return;
511   }
512}
513
514double log(const xdouble& a)
515{
516   static double LogBound = log(NTL_XD_BOUND);
517   if (a.x <= 0) {
518      Error("log(xdouble): argument must be positive");
519   }
520
521   return log(a.x) + a.e*LogBound;
522}
523
524xdouble xexp(double x)
525{
526   const double LogBound = log(NTL_XD_BOUND);
527
528   double y = x/LogBound;
529   double iy = floor(y+0.5);
530
531   if (iy >= NTL_OVFBND)
532      Error("xdouble: overflow");
533
534   if (iy <= -NTL_OVFBND)
535      Error("xdouble: underflow");
536
537
538   double fy = y - iy;
539
540   xdouble res;
541   res.e = long(iy);
542   res.x = exp(fy*LogBound);
543   res.normalize();
544   return res;
545}
546
547/**************  input / output routines **************/
548
549
550void ComputeLn2(RR&);
551void ComputeLn10(RR&);
552
553long ComputeMax10Power()
554{
555   long old_p = RR::precision();
556   RR::SetPrecision(NTL_BITS_PER_LONG);
557
558   RR ln2, ln10;
559   ComputeLn2(ln2);
560   ComputeLn10(ln10);
561
562   long k = to_long( to_RR(NTL_OVFBND/2) * ln2 / ln10 );
563
564   RR::SetPrecision(old_p);
565   return k;
566}
567
568
569xdouble PowerOf10(const ZZ& e)
570{
571   static long init = 0;
572   static xdouble v10k;
573   static long k;
574
575   if (!init) {
576      long old_p = RR::precision();
577      k = ComputeMax10Power();
578      RR::SetPrecision(NTL_DOUBLE_PRECISION);
579      v10k = to_xdouble(power(to_RR(10), k)); 
580      RR::SetPrecision(old_p);
581      init = 1;
582   }
583
584   ZZ e1;
585   long neg;
586
587   if (e < 0) {
588      e1 = -e;
589      neg = 1;
590   }
591   else {
592      e1 = e;
593      neg = 0;
594   }
595
596   long r;
597   ZZ q;
598
599   r = DivRem(q, e1, k);
600
601   long old_p = RR::precision();
602   RR::SetPrecision(NTL_DOUBLE_PRECISION);
603   xdouble x1 = to_xdouble(power(to_RR(10), r));
604   RR::SetPrecision(old_p);
605
606   xdouble x2 = power(v10k, q);
607   xdouble x3 = x1*x2;
608
609   if (neg) x3 = 1/x3;
610
611   return x3;
612}
613
614
615
616
617xdouble to_xdouble(const char *s)
618{
619   long c;
620   long cval;
621   long sign;
622   ZZ a, b;
623   long i=0;
624
625   if (!s) Error("bad xdouble input");
626
627   c = s[i];
628   while (IsWhiteSpace(c)) {
629      i++;
630      c = s[i];
631   }
632
633   if (c == '-') {
634      sign = -1;
635      i++;
636      c = s[i];
637   }
638   else
639      sign = 1;
640
641   long got1 = 0;
642   long got_dot = 0;
643   long got2 = 0;
644
645   a = 0;
646   b = 1;
647
648   cval = CharToIntVal(c);
649
650   if (cval >= 0 && cval <= 9) {
651      got1 = 1;
652
653      while (cval >= 0 && cval <= 9) {
654         mul(a, a, 10);
655         add(a, a, cval);
656         i++;
657         c = s[i];
658         cval = CharToIntVal(c);
659      }
660   }
661
662   if (c == '.') {
663      got_dot = 1;
664
665      i++;
666      c = s[i];
667      cval = CharToIntVal(c);
668
669      if (cval >= 0 && cval <= 9) {
670         got2 = 1;
671   
672         while (cval >= 0 && cval <= 9) {
673            mul(a, a, 10);
674            add(a, a, cval);
675            mul(b, b, 10);
676            i++;
677            c = s[i];
678            cval = CharToIntVal(c);
679         }
680      }
681   }
682
683   if (got_dot && !got1 && !got2)  Error("bad xdouble input");
684
685   ZZ e;
686
687   long got_e = 0;
688   long e_sign;
689
690   if (c == 'e' || c == 'E') {
691      got_e = 1;
692
693      i++;
694      c = s[i];
695
696      if (c == '-') {
697         e_sign = -1;
698         i++;
699         c = s[i];
700      }
701      else if (c == '+') {
702         e_sign = 1;
703         i++;
704         c = s[i];
705      }
706      else
707         e_sign = 1;
708
709      cval = CharToIntVal(c);
710
711      if (cval < 0 || cval > 9) Error("bad xdouble input");
712
713      e = 0;
714      while (cval >= 0 && cval <= 9) {
715         mul(e, e, 10);
716         add(e, e, cval);
717         i++;
718         c = s[i];
719         cval = CharToIntVal(c);
720      }
721   }
722
723   if (!got1 && !got2 && !got_e) Error("bad xdouble input");
724
725   xdouble t1, t2, v;
726
727   if (got1 || got2) {
728      conv(t1, a);
729      conv(t2, b);
730      v = t1/t2;
731   }
732   else
733      v = 1;
734
735   if (sign < 0)
736      v = -v;
737
738   if (got_e) {
739      if (e_sign < 0) negate(e, e);
740      t1 = PowerOf10(e);
741      v = v * t1;
742   }
743
744   return v;
745}
746
747NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.