source: git/ntl/src/RR.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: 27.0 KB
Line 
1
2#include <NTL/RR.h>
3
4
5#include <NTL/new.h>
6
7NTL_START_IMPL
8
9
10long RR::prec = 150;
11
12void RR::SetPrecision(long p)
13{
14   if (p < 53)
15      p = 53;
16
17   if (NTL_OVERFLOW(p, 1, 0))
18      Error("RR: precision too high");
19
20   prec = p;
21}
22
23long RR::oprec = 10;
24
25void RR::SetOutputPrecision(long p)
26{
27   if (p < 1)
28      p = 1;
29
30   if (NTL_OVERFLOW(p, 1, 0))
31      Error("RR: output precision too high");
32
33   oprec = p;
34}
35
36
37
38static
39void normalize1(RR& z, const ZZ& y_x, long y_e, long prec, long residual)
40{
41   long len = NumBits(y_x);
42
43   if (len > prec) {
44      long correction = ZZ_RoundCorrection(y_x, len - prec, residual);
45
46      RightShift(z.x, y_x, len - prec);
47
48      if (correction) 
49         add(z.x, z.x, correction);
50
51      z.e = y_e + len - prec;
52   }
53   else if (len == 0) {
54      clear(z.x);
55      z.e = 0;
56   }
57   else {
58      z.x = y_x;
59      z.e = y_e;
60   }
61
62   if (!IsOdd(z.x))
63      z.e += MakeOdd(z.x);
64
65   if (z.e >= NTL_OVFBND)
66      Error("RR: overflow");
67
68   if (z.e <= -NTL_OVFBND)
69      Error("RR: underflow");
70}
71
72void normalize(RR& z, const RR& y, long residual = 0)
73{
74   normalize1(z, y.x, y.e, RR::prec, residual);
75}
76
77void MakeRR(RR& z, const ZZ& a,  long e)
78{
79   if (e >= NTL_OVFBND)
80      Error("MakeRR: e too big");
81
82   if (e <= -NTL_OVFBND)
83      Error("MakeRR: e too small");
84
85   normalize1(z, a, e, RR::prec, 0);
86}
87
88void MakeRRPrec(RR& x, const ZZ& a, long e, long p)
89{
90   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
91      Error("MakeRRPrec: bad precsion");
92
93   long old_p = RR::prec;
94   RR::prec = p;
95   MakeRR(x, a, e);
96   RR::prec = old_p;
97}
98
99void random(RR& z)
100{
101   static RR t;
102   RandomBits(t.x, RR::prec); 
103   t.e = -RR::prec;
104   normalize(z, t);
105}
106
107
108static inline 
109void xcopy(RR& x, const RR& a)
110   { normalize(x, a); }
111
112// xcopy emulates old assignment semantics...
113// many routines here implicitly assume assignment normalizes,
114// but that is no longer the case as of v3.0.
115
116void ConvPrec(RR& x, const RR& a, long p)
117{
118   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
119      Error("ConvPrec: bad precsion");
120
121   long old_p = RR::prec;
122   RR::prec = p;
123   normalize(x, a);
124   RR::prec = old_p;
125}
126
127void RoundToPrecision(RR& x, const RR& a, long p)
128{
129   ConvPrec(x, a, p);
130}
131   
132
133void conv(RR& x, const RR& a)
134{
135   normalize(x, a);
136}
137
138
139long IsZero(const RR& a)
140{
141   return IsZero(a.x);
142}
143
144long IsOne(const RR& a)
145{
146   return a.e == 0 && IsOne(a.x);
147}
148
149long sign(const RR& a)
150{
151   return sign(a.x);
152}
153
154void clear(RR& z)
155{
156   z.e = 0;
157   clear(z.x);
158}
159
160void set(RR& z)
161{
162   z.e = 0;
163   set(z.x);
164}
165
166
167void add(RR& z, const RR& a, const RR& b)
168{
169   static RR t;
170
171   if (IsZero(a.x)) {
172      xcopy(z, b);
173      return;
174   }
175
176   if (IsZero(b.x)) {
177      xcopy(z, a);
178      return;
179   }
180
181   if (a.e > b.e) {
182      if (a.e-b.e - max(RR::prec-NumBits(a.x),0) >= NumBits(b.x) + 2)
183         normalize(z, a, sign(b));
184      else {
185         LeftShift(t.x, a.x, a.e-b.e);
186         add(t.x, t.x, b.x);
187         t.e = b.e;
188         normalize(z, t);
189      }
190   }
191   else if (a.e < b.e) {
192      if (b.e-a.e - max(RR::prec-NumBits(b.x),0) >= NumBits(a.x) + 2)
193         normalize(z, b, sign(a));
194      else {
195         LeftShift(t.x, b.x, b.e-a.e);
196         add(t.x, t.x, a.x);
197         t.e = a.e;
198         normalize(z, t);
199      }
200   }
201   else {
202      add(t.x, a.x, b.x);
203      t.e = a.e;
204      normalize(z, t);
205   }
206}
207
208void AddPrec(RR& x, const RR& a, const RR& b, long p)
209{
210   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
211      Error("AddPrec: bad precsion");
212
213   long old_p = RR::prec;
214   RR::prec = p;
215   add(x, a, b);
216   RR::prec = old_p;
217}
218
219void sub(RR& z, const RR& a, const RR& b)
220{
221   static RR t;
222
223   if (IsZero(a.x)) {
224      negate(z, b);
225      return;
226   }
227
228   if (IsZero(b.x)) {
229      xcopy(z, a);
230      return;
231   }
232
233   if (a.e > b.e) {
234      if (a.e-b.e - max(RR::prec-NumBits(a.x),0) >= NumBits(b.x) + 2)
235         normalize(z, a, -sign(b));
236      else {
237         LeftShift(t.x, a.x, a.e-b.e);
238         sub(t.x, t.x, b.x);
239         t.e = b.e;
240         xcopy(z, t);
241      }
242   }
243   else if (a.e < b.e) {
244      if (b.e-a.e - max(RR::prec-NumBits(b.x),0) >= NumBits(a.x) + 2) {
245         normalize(z, b, -sign(a));
246         negate(z.x, z.x);
247      }
248      else {
249         LeftShift(t.x, b.x, b.e-a.e);
250         sub(t.x, a.x, t.x);
251         t.e = a.e;
252         xcopy(z, t);
253      }
254   }
255   else {
256      sub(t.x, a.x, b.x);
257      t.e = a.e;
258      normalize(z, t);
259   }
260}
261
262void SubPrec(RR& x, const RR& a, const RR& b, long p)
263{
264   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
265      Error("SubPrec: bad precsion");
266
267   long old_p = RR::prec;
268   RR::prec = p;
269   sub(x, a, b);
270   RR::prec = old_p;
271}
272
273void negate(RR& z, const RR& a)
274{
275   xcopy(z, a);
276   negate(z.x, z.x);
277}
278
279void NegatePrec(RR& x, const RR& a, const RR& b, long p)
280{
281   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
282      Error("NegatePrec: bad precsion");
283
284   long old_p = RR::prec;
285   RR::prec = p;
286   negate(x, a);
287   RR::prec = old_p;
288}
289
290void abs(RR& z, const RR& a)
291{
292   xcopy(z, a);
293   abs(z.x, z.x);
294}
295
296void AbsPrec(RR& x, const RR& a, const RR& b, long p)
297{
298   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
299      Error("AbsPrec: bad precsion");
300
301   long old_p = RR::prec;
302   RR::prec = p;
303   abs(x, a);
304   RR::prec = old_p;
305}
306
307
308void mul(RR& z, const RR& a, const RR& b)
309{
310   static RR t;
311
312   mul(t.x, a.x, b.x);
313   t.e = a.e + b.e;
314   xcopy(z, t);
315}
316
317void MulPrec(RR& x, const RR& a, const RR& b, long p)
318{
319   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
320      Error("MulPrec: bad precsion");
321
322   long old_p = RR::prec;
323   RR::prec = p;
324   mul(x, a, b);
325   RR::prec = old_p;
326}
327
328
329void sqr(RR& z, const RR& a)
330{
331   static RR t;
332
333   sqr(t.x, a.x);
334   t.e = a.e + a.e;
335   xcopy(z, t);
336}
337
338void SqrPrec(RR& x, const RR& a, const RR& b, long p)
339{
340   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
341      Error("SqrPrec: bad precsion");
342
343   long old_p = RR::prec;
344   RR::prec = p;
345   sqr(x, a);
346   RR::prec = old_p;
347}
348
349
350
351void div(RR& z, const RR& a, const RR& b)
352{
353   if (IsZero(b))
354      Error("RR: division by zero");
355
356   if (IsZero(a)) {
357      clear(z);
358      return;
359   }
360
361   long la = NumBits(a.x);
362   long lb = NumBits(b.x);
363
364   long neg = (sign(a) != sign(b));
365
366   long k = RR::prec - la + lb + 1;
367   if (k < 0) k = 0;
368
369   static RR t;
370   static ZZ A, B, R;
371
372   abs(A, a.x);
373   LeftShift(A, A, k);
374
375   abs(B, b.x);
376   DivRem(t.x, R, A, B);
377
378   t.e = a.e - b.e - k;
379
380   normalize(z, t, !IsZero(R));
381
382   if (neg)
383      negate(z.x, z.x);
384}
385
386void DivPrec(RR& x, const RR& a, const RR& b, long p)
387{
388   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
389      Error("DivPrec: bad precsion");
390
391   long old_p = RR::prec;
392   RR::prec = p;
393   div(x, a, b);
394   RR::prec = old_p;
395}
396
397void SqrRoot(RR& z, const RR& a)
398{
399   if (sign(a) < 0)
400      Error("RR: attempt to take square root of negative number");
401
402   if (IsZero(a)) {
403      clear(z);
404      return;
405   }
406
407   RR t;
408   ZZ T1, T2;
409   long k;
410
411   k = 2*RR::prec - NumBits(a.x) + 1;
412
413   if (k < 0) k = 0;
414
415   if ((a.e - k) & 1) k++;
416
417   LeftShift(T1, a.x, k);
418   SqrRoot(t.x, T1);
419   t.e = (a.e - k)/2;
420   sqr(T2, T1);
421
422   normalize(z, t, T2 < T1);
423}
424
425void SqrRootPrec(RR& x, const RR& a, const RR& b, long p)
426{
427   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
428      Error("SqrRootPrec: bad precsion");
429
430   long old_p = RR::prec;
431   RR::prec = p;
432   SqrRoot(x, a);
433   RR::prec = old_p;
434}
435
436
437
438
439void swap(RR& a, RR& b)
440{
441   swap(a.x, b.x);
442   swap(a.e, b.e);
443}
444
445long compare(const RR& a, const RR& b)
446{
447   static RR t;
448
449   SubPrec(t, a, b, 1);
450   return sign(t);
451}
452
453
454
455long operator==(const RR& a, const RR& b) 
456{
457   return a.e == b.e && a.x == b.x;
458}
459
460
461void trunc(RR& z, const RR& a)
462{
463   static RR t;
464
465   if (a.e >= 0) 
466      xcopy(z, a);
467   else {
468      RightShift(t.x, a.x, -a.e);
469      t.e = 0;
470      xcopy(z, t);
471   }
472}
473
474void TruncPrec(RR& x, const RR& a, const RR& b, long p)
475{
476   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
477      Error("TruncPrec: bad precsion");
478
479   long old_p = RR::prec;
480   RR::prec = p;
481   trunc(x, a);
482   RR::prec = old_p;
483}
484
485void floor(RR& z, const RR& a)
486{
487   static RR t;
488
489   if (a.e >= 0) 
490      xcopy(z, a);
491   else {
492      RightShift(t.x, a.x, -a.e);
493      if (sign(a.x) < 0)
494         add(t.x, t.x, -1);
495      t.e = 0;
496      xcopy(z, t);
497   }
498}
499
500void FloorPrec(RR& x, const RR& a, const RR& b, long p)
501{
502   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
503      Error("FloorPrec: bad precsion");
504
505   long old_p = RR::prec;
506   RR::prec = p;
507   floor(x, a);
508   RR::prec = old_p;
509}
510
511void ceil(RR& z, const RR& a)
512{
513   static RR t;
514
515   if (a.e >= 0)
516      xcopy(z, a);
517   else {
518      RightShift(t.x, a.x, -a.e);
519      if (sign(a.x) > 0)
520         add(t.x, t.x, 1);
521      t.e = 0;
522      xcopy(z, t);
523   }
524}
525
526void CeilPrec(RR& x, const RR& a, const RR& b, long p)
527{
528   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
529      Error("CeilPrec: bad precsion");
530
531   long old_p = RR::prec;
532   RR::prec = p;
533   ceil(x, a);
534   RR::prec = old_p;
535}
536
537void round(RR& z, const RR& a)
538{
539   if (a.e >= 0) {
540      xcopy(z, a);
541      return;
542   }
543
544   long len = NumBits(a.x);
545
546   if (-a.e > len) {
547      z = 0;
548      return;
549   }
550
551   if (-a.e == len) {
552      if (len == 1)
553         z = 0;
554      else
555         z = sign(a.x);
556
557      return;
558   }
559
560   static RR t;
561   ConvPrec(t, a, len+a.e);
562   xcopy(z, t);
563}
564
565void RoundPrec(RR& x, const RR& a, const RR& b, long p)
566{
567   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
568      Error("RoundPrec: bad precsion");
569
570   long old_p = RR::prec;
571   RR::prec = p;
572   round(x, a);
573   RR::prec = old_p;
574}
575
576
577   
578
579void conv(RR& z, const ZZ& a)
580{
581   normalize1(z, a, 0, RR::prec, 0);
582}
583
584void ConvPrec(RR& x, const ZZ& a, long p)
585{
586   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
587      Error("ConvPrec: bad precsion");
588
589   long old_p = RR::prec;
590   RR::prec = p;
591   conv(x, a);
592   RR::prec = old_p;
593}
594
595
596void conv(RR& z, long a)
597{
598   if (a == 0) {
599      clear(z);
600      return;
601   }
602
603   if (a == 1) {
604      set(z);
605      return;
606   }
607
608   static ZZ t;
609   t = a;
610   conv(z, t);
611}
612
613void ConvPrec(RR& x, long a, long p)
614{
615   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
616      Error("ConvPrec: bad precsion");
617
618   long old_p = RR::prec;
619   RR::prec = p;
620   conv(x, a);
621   RR::prec = old_p;
622}
623
624void conv(RR& z, unsigned long a)
625{
626   if (a == 0) {
627      clear(z);
628      return;
629   }
630
631   if (a == 1) {
632      set(z);
633      return;
634   }
635
636   static ZZ t;
637   conv(t, a);
638   conv(z, t);
639}
640
641void ConvPrec(RR& x, unsigned long a, long p)
642{
643   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
644      Error("ConvPrec: bad precsion");
645
646   long old_p = RR::prec;
647   RR::prec = p;
648   conv(x, a);
649   RR::prec = old_p;
650}
651
652
653void conv(RR& z, double a)
654{
655   if (a == 0) {
656      clear(z);
657      return;
658   }
659
660   if (a == 1) {
661      set(z);
662      return;
663   }
664
665   if (!IsFinite(&a))
666      Error("RR: conversion of a non-finite double");
667
668   int e;
669   double f;
670   static RR t;
671
672   f = frexp(a, &e);
673
674   f = f * NTL_FDOUBLE_PRECISION;
675   f = f * 4;
676
677   conv(t.x, f);
678   t.e = e - (NTL_DOUBLE_PRECISION + 1);
679
680   xcopy(z, t);
681}
682
683void ConvPrec(RR& x, double a, long p)
684{
685   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
686      Error("ConvPrec: bad precsion");
687
688   long old_p = RR::prec;
689   RR::prec = p;
690   conv(x, a);
691   RR::prec = old_p;
692}
693
694
695void conv(ZZ& z, const RR& a)
696{
697   if (a.e >= 0) 
698      LeftShift(z, a.x, a.e);
699   else {
700      long sgn = sign(a.x);
701      RightShift(z, a.x, -a.e);
702      if (sgn < 0)
703         sub(z, z, 1);
704   }
705}
706
707void CeilToZZ(ZZ& z, const RR& a)
708{
709   if (a.e >= 0)
710      LeftShift(z, a.x, a.e);
711   else {
712      long sgn = sign(a.x);
713      RightShift(z, a.x, -a.e);
714      if (sgn > 0)
715         add(z, z, 1);
716   }
717}
718
719void TruncToZZ(ZZ& z, const RR& a)
720{
721   if (a.e >= 0)
722      LeftShift(z, a.x, a.e);
723   else 
724      RightShift(z, a.x, -a.e);
725}
726
727
728void RoundToZZ(ZZ& z, const RR& a)
729{
730   if (a.e >= 0) {
731      LeftShift(z, a.x, a.e);
732      return;
733   }
734
735   long len = NumBits(a.x);
736
737   if (-a.e > len) {
738      z = 0;
739      return;
740   }
741
742   if (-a.e == len) {
743      if (len == 1)
744         z = 0;
745      else
746         z = sign(a.x);
747
748      return;
749   }
750
751   static RR t;
752
753   ConvPrec(t, a, len+a.e);
754
755   LeftShift(z, t.x, t.e);
756}
757
758
759void conv(long& z, const RR& a)
760{
761   ZZ t;
762   conv(t, a);
763   conv(z, t);
764}
765
766void conv(double& z, const RR& aa)
767{
768   double x;
769   static RR a;
770
771   ConvPrec(a, aa, NTL_DOUBLE_PRECISION);
772   // round to NTL_DOUBLE_PRECISION bits to avoid double overflow
773
774   conv(x, a.x);
775   z = _ntl_ldexp(x, a.e);
776}
777
778
779void add(RR& z, const RR& a, double b)
780{
781   static RR B;
782   B = b;
783   add(z, a, B);
784}
785
786
787
788void sub(RR& z, const RR& a, double b)
789{
790   static RR B;
791   B = b;
792   sub(z, a, B);
793}
794
795void sub(RR& z, double a, const RR& b)
796{
797   static RR A;
798   A = a;
799   sub(z, A, b);
800}
801
802
803
804void mul(RR& z, const RR& a, double b)
805{
806   static RR B;
807   B = b;
808   mul(z, a, B);
809}
810
811
812void div(RR& z, const RR& a, double b)
813{
814   static RR B;
815   B = b;
816   div(z, a, B);
817}
818
819void div(RR& z, double a, const RR& b)
820{
821   static RR A;
822   A = a;
823   div(z, A, b);
824}
825
826
827void inv(RR& z, const RR& a)
828{
829   static RR one = to_RR(1);
830   div(z, one, a);
831}
832
833void InvPrec(RR& x, const RR& a, long p)
834{
835   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
836      Error("InvPrec: bad precsion");
837
838   long old_p = RR::prec;
839   RR::prec = p;
840   inv(x, a);
841   RR::prec = old_p;
842}
843
844
845long compare(const RR& a, double b)
846{
847   if (b == 0) return sign(a);
848
849   static RR B;
850   B = b;
851   return compare(a, B);
852}
853
854
855long operator==(const RR& a, double b) 
856{
857   if (b == 0) return IsZero(a);
858   if (b == 1) return IsOne(a);
859
860   static RR B;
861   B = b;
862   return a == B;
863}
864
865
866void power(RR& z, const RR& a, long e)
867{
868   RR b, res;
869
870   long n = NumBits(e);
871
872   long p = RR::precision();
873   RR::SetPrecision(p + n + 10);
874
875   xcopy(b, a);
876
877   set(res);
878   long i;
879
880   for (i = n-1; i >= 0; i--) {
881      sqr(res, res);
882      if (bit(e, i))
883         mul(res, res, b);
884   }
885
886   RR::SetPrecision(p);
887
888   if (e < 0) 
889      inv(z, res);
890   else
891      xcopy(z, res);
892}
893
894
895void conv(RR& z, const xdouble& a)
896{
897   conv(z, a.mantissa());
898
899   if (a.exponent() >  ((2*NTL_OVFBND)/(2*NTL_XD_HBOUND_LOG))) 
900      Error("RR: overlow");
901
902   if (a.exponent() < -((2*NTL_OVFBND)/(2*NTL_XD_HBOUND_LOG))) 
903      Error("RR: underflow");
904
905   z.e += a.exponent()*(2*NTL_XD_HBOUND_LOG);
906
907   if (z.e >= NTL_OVFBND)
908      Error("RR: overflow");
909
910   if (z.e <= -NTL_OVFBND)
911      Error("RR: underflow");
912}
913
914void ConvPrec(RR& x, const xdouble& a, long p)
915{
916   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
917      Error("ConvPrec: bad precsion");
918
919   long old_p = RR::prec;
920   RR::prec = p;
921   conv(x, a);
922   RR::prec = old_p;
923}
924
925
926
927void conv(xdouble& z, const RR& a)
928{
929   xdouble x;
930   xdouble y;
931
932   conv(x, a.x);
933   power2(y, a.e);
934   z = x*y;
935}
936     
937void power2(RR& z, long e)
938{
939   if (e >= NTL_OVFBND)
940      Error("RR: overflow");
941
942   if (e <= -NTL_OVFBND)
943      Error("RR: underflow");
944
945   set(z.x); 
946   z.e = e;
947}
948
949void conv(RR& z, const quad_float& a)
950{
951   static RR hi, lo, res;
952
953   ConvPrec(hi, a.hi, NTL_DOUBLE_PRECISION);
954   ConvPrec(lo, a.lo, NTL_DOUBLE_PRECISION);
955
956   add(res, hi, lo);
957
958   z = res;
959}
960
961void ConvPrec(RR& x, const quad_float& a, long p)
962{
963   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
964      Error("ConvPrec: bad precsion");
965
966   long old_p = RR::prec;
967   RR::prec = p;
968   conv(x, a);
969   RR::prec = old_p;
970}
971
972
973void conv(quad_float& z, const RR& a)
974{
975   long old_p;
976   static RR a_hi, a_lo;
977
978   old_p = RR::prec;
979
980   ConvPrec(a_hi, a, NTL_DOUBLE_PRECISION);  // high order bits
981   SubPrec(a_lo, a, a_hi, NTL_DOUBLE_PRECISION);  // low order bits
982
983   z = to_quad_float(a_hi.x)*power2_quad_float(a_hi.e) +
984       to_quad_float(a_lo.x)*power2_quad_float(a_lo.e);
985}
986
987void conv(RR& x, const char *s)
988{
989   long c;
990   long cval;
991   long sign;
992   ZZ a, b;
993   long i = 0;
994
995   if (!s) Error("bad RR input");
996
997
998   c = s[i];
999   while (IsWhiteSpace(c)) {
1000      i++;
1001      c = s[i];
1002   }
1003
1004   if (c == '-') {
1005      sign = -1;
1006      i++;
1007      c = s[i];
1008   }
1009   else
1010      sign = 1;
1011
1012   long got1 = 0;
1013   long got_dot = 0;
1014   long got2 = 0;
1015
1016   a = 0;
1017   b = 1;
1018
1019   cval = CharToIntVal(c);
1020
1021   if (cval >= 0 && cval <= 9) {
1022      got1 = 1;
1023
1024      while (cval >= 0 && cval <= 9) {
1025         mul(a, a, 10);
1026         add(a, a, cval);
1027         i++;
1028         c = s[i];
1029         cval = CharToIntVal(c);
1030      }
1031   }
1032
1033   if (c == '.') {
1034      got_dot = 1;
1035
1036      i++;
1037      c = s[i];
1038      cval = CharToIntVal(c);
1039
1040      if (cval >= 0 && cval <= 9) {
1041         got2 = 1;
1042   
1043         while (cval >= 0 && cval <= 9) {
1044            mul(a, a, 10);
1045            add(a, a, cval);
1046            mul(b, b, 10);
1047            i++;
1048            c = s[i];
1049            cval = CharToIntVal(c);
1050         }
1051      }
1052   }
1053
1054   if (got_dot && !got1 && !got2)  Error("bad RR input");
1055
1056   ZZ e;
1057
1058   long got_e = 0;
1059   long e_sign;
1060
1061   if (c == 'e' || c == 'E') {
1062      got_e = 1;
1063
1064      i++;
1065      c = s[i];
1066
1067      if (c == '-') {
1068         e_sign = -1;
1069         i++;
1070         c = s[i];
1071      }
1072      else if (c == '+') {
1073         e_sign = 1;
1074         i++;
1075         c = s[i];
1076      }
1077      else
1078         e_sign = 1;
1079
1080
1081      cval = CharToIntVal(c);
1082
1083      if (cval < 0 || cval > 9) Error("bad RR input");
1084
1085      e = 0;
1086      while (cval >= 0 && cval <= 9) {
1087         mul(e, e, 10);
1088         add(e, e, cval);
1089         i++;
1090         c = s[i];
1091         cval = CharToIntVal(c);
1092      }
1093   }
1094
1095   if (!got1 && !got2 && !got_e) Error("bad RR input");
1096
1097   RR t1, t2, v;
1098
1099   long old_p = RR::precision();
1100
1101   if (got1 || got2) {
1102      ConvPrec(t1, a, max(NumBits(a), 1));
1103      ConvPrec(t2, b, NumBits(b));
1104      if (got_e)
1105         RR::SetPrecision(old_p + 10);
1106
1107      div(v, t1, t2);
1108   }
1109   else
1110      set(v);
1111
1112   if (sign < 0)
1113      negate(v, v);
1114
1115   if (got_e) {
1116      if (e >= NTL_OVFBND) Error("RR input overflow");
1117      long E;
1118      conv(E, e);
1119      if (e_sign < 0) E = -E;
1120      RR::SetPrecision(old_p + 10);
1121      power(t1, to_RR(10), E);
1122      mul(v, v, t1);
1123      RR::prec = old_p;
1124   }
1125
1126   xcopy(x, v);
1127}
1128
1129void ConvPrec(RR& x, const char *s, long p)
1130{
1131   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
1132      Error("ConvPrec: bad precsion");
1133
1134   long old_p = RR::prec;
1135   RR::prec = p;
1136   conv(x, s);
1137   RR::prec = old_p;
1138}
1139
1140
1141void ReallyComputeE(RR& res)
1142{
1143   long p = RR::precision();
1144   RR::SetPrecision(p + NumBits(p) + 10);
1145
1146   RR s, s1, t;
1147
1148   s = 1;
1149   t = 1;
1150
1151   long i;
1152
1153   for (i = 2; ; i++) {
1154      add(s1, s, t);
1155      if (s == s1) break;
1156      xcopy(s, s1);
1157      div(t, t, i);
1158   }
1159
1160   RR::SetPrecision(p);
1161   xcopy(res, s);
1162}
1163
1164void ComputeE(RR& res)
1165{
1166   static long prec = 0;
1167   static RR e;
1168
1169   long p = RR::precision();
1170
1171   if (prec <= p + 10) {
1172      prec = p + 20;
1173      RR::SetPrecision(prec);
1174      ReallyComputeE(e);
1175      RR::SetPrecision(p);
1176   }
1177
1178   xcopy(res, e);
1179}
1180
1181
1182void exp(RR& res, const RR& x)
1183{
1184   if (x >= NTL_OVFBND || x <= -NTL_OVFBND)
1185      Error("RR: overflow");
1186
1187   long p = RR::precision();
1188
1189   // step 0: write x = n + f, n an integer and |f| <= 1/2
1190   //    careful -- we want to compute f to > p bits of precision
1191
1192
1193   RR f, nn;
1194   RR::SetPrecision(NTL_BITS_PER_LONG);
1195   round(nn, x);
1196   RR::SetPrecision(p + 10);
1197   sub(f, x, nn);
1198   long n = to_long(nn);
1199
1200   // step 1: calculate t1 = e^n by repeated squaring
1201
1202   RR::SetPrecision(p + NumBits(n) + 10);
1203
1204   RR e;
1205   ComputeE(e);
1206
1207   RR::SetPrecision(p + 10);
1208
1209   RR t1;
1210   power(t1, e, n); 
1211
1212   // step 2: calculate t2 = e^f using Taylor series expansion
1213
1214   RR::SetPrecision(p + NumBits(p) + 10);
1215
1216   RR t2, s, s1, t;
1217   long i;
1218
1219   s = 0;
1220   t = 1;
1221
1222   for (i = 1; ; i++) {
1223      add(s1, s, t);
1224      if (s == s1) break;
1225      xcopy(s, s1);
1226      mul(t, t, f);
1227      div(t, t, i);
1228   }
1229
1230   xcopy(t2, s);
1231
1232   RR::SetPrecision(p);
1233
1234   mul(res, t1, t2);
1235}
1236
1237
1238
1239void ReallyComputeLn2(RR& res)
1240{
1241   long p = RR::precision();
1242   RR::SetPrecision(p + NumBits(p) + 10);
1243
1244   RR s, s1, t, t1;
1245
1246   s = 0;
1247   t = 0.5;
1248   t1 = 0.5;
1249
1250   long i;
1251
1252   for (i = 2; ; i++) {
1253      add(s1, s, t);
1254      if (s == s1) break;
1255      xcopy(s, s1);
1256      mul(t1, t1, 0.5);
1257      div(t, t1, i);
1258   }
1259
1260   RR::SetPrecision(p);
1261   xcopy(res, s);
1262}
1263
1264
1265void ComputeLn2(RR& res)
1266{
1267   static long prec = 0;
1268   static RR ln2;
1269
1270   long p = RR::precision();
1271
1272   if (prec <= p + 10) {
1273      prec = p + 20;
1274      RR::SetPrecision(prec);
1275      ReallyComputeLn2(ln2);
1276      RR::SetPrecision(p);
1277   }
1278
1279   xcopy(res, ln2);
1280}
1281
1282long Lg2(const RR& x)
1283{
1284   return NumBits(x.mantissa()) + x.exponent();
1285}
1286
1287void log(RR& res, const RR& x)
1288{
1289   if (x <= 0) Error("argument to log must be positive");
1290
1291   long p = RR::precision();
1292
1293   RR::SetPrecision(p + NumBits(p) + 10);
1294
1295   RR y;
1296   long n;
1297
1298   // re-write x = 2^n * (1 - y), where -1/2 < y < 1/4  (so 3/4 < 1-y < 3/2)
1299
1300   if (x > 0.75 && x < 1.5) {
1301      n = 0;
1302      sub(y, 1, x);
1303   }
1304   else {
1305      n = Lg2(x) - 1;
1306      RR t;
1307      power2(t, -n);
1308      mul(t, t, x);
1309      while (t > 1.5) {
1310         mul(t, t, 0.5);
1311         n++;
1312      }
1313
1314      sub(y, 1, t);
1315   }
1316
1317   // compute s = - ln(1-y) by power series expansion
1318
1319   RR s, s1, t, t1;
1320
1321   s = 0;
1322   xcopy(t, y);
1323   xcopy(t1, y);
1324
1325   long i;
1326
1327   for (i = 2; ; i++) {
1328      add(s1, s, t);
1329      if (s == s1) break;
1330      xcopy(s, s1);
1331      mul(t1, t1, y);
1332      div(t, t1, i);
1333   }
1334
1335   if (n == 0) 
1336      t = 0;
1337   else {
1338      ComputeLn2(t);
1339      mul(t, t, n);
1340   }
1341
1342   RR::SetPrecision(p);
1343
1344   sub(res, t, s);
1345}
1346
1347
1348void ComputeLn10(RR& res)
1349{
1350   static long prec = 0;
1351   static RR ln10;
1352
1353   long p = RR::precision();
1354
1355   if (prec <= p + 10) {
1356      prec = p + 20;
1357      RR::SetPrecision(prec);
1358      log(ln10, to_RR(10));
1359      RR::SetPrecision(p);
1360   }
1361
1362   xcopy(res, ln10);
1363}
1364
1365void log10(RR& res, const RR& x)
1366{
1367   long p = RR::precision();
1368   RR::SetPrecision(p + 10);
1369
1370   RR ln10, t1, t2;
1371   ComputeLn10(ln10);
1372
1373   log(t1, x);
1374   div(t2, t1, ln10);
1375
1376   RR::SetPrecision(p);
1377
1378   xcopy(res, t2);
1379}
1380
1381
1382void expm1(RR& res, const RR& x)
1383{
1384   if (x < -0.5 || x > 0.5) {
1385      RR t;
1386      exp(t, x);
1387      sub(res, t, 1);
1388      return;
1389   }
1390
1391   long p = RR::precision();
1392
1393   RR::SetPrecision(p + NumBits(p) + 10);
1394
1395   RR f;
1396
1397   xcopy(f, x);
1398
1399   RR s, s1, t;
1400   long i;
1401
1402   s = 0;
1403   xcopy(t, f);
1404
1405   for (i = 2; ; i++) {
1406      add(s1, s, t);
1407      if (s == s1) break;
1408      xcopy(s, s1);
1409      mul(t, t, f);
1410      div(t, t, i);
1411   }
1412
1413   RR::SetPrecision(p);
1414
1415   xcopy(res, s);
1416}
1417
1418
1419
1420void log1p(RR& res, const RR& x)
1421{
1422   if (x < -0.5 || x > 0.5) {
1423      log(res, 1 + x);
1424      return;
1425   }
1426
1427   long p = RR::precision();
1428
1429   RR::SetPrecision(p + NumBits(p) + 10);
1430
1431   RR y;
1432
1433   negate(y, x);
1434
1435   // compute s = - ln(1-y) by power series expansion
1436
1437   RR s, s1, t, t1;
1438
1439   s = 0;
1440   xcopy(t, y);
1441   xcopy(t1, y);
1442
1443   long i;
1444
1445   for (i = 2; ; i++) {
1446      add(s1, s, t);
1447      if (s == s1) break;
1448      xcopy(s, s1);
1449      mul(t1, t1, y);
1450      div(t, t1, i);
1451   }
1452
1453   RR::SetPrecision(p);
1454
1455   negate(res, s);
1456
1457}
1458
1459
1460void pow(RR& res, const RR& x, const RR& y)
1461{
1462   if (y == 0) {
1463      res = 1;
1464      return;
1465   }
1466
1467   if (x == 0) {
1468      res = 0;
1469      return;
1470   }
1471
1472   if (x == 1) {
1473      res = 1;
1474      return;
1475   }
1476
1477   if (x < 0) {
1478      Error("pow: sorry...first argument to pow must be nonnegative");
1479   }
1480
1481   long p = RR::precision();
1482
1483   // calculate working precison...one could use p + NTL_BITS_PER_LONG + 10,
1484   // for example, but we want the behaviour to be machine independent.
1485   // so we calculate instead a rough approximation to log |y log(x)|
1486
1487   RR t1, t2;
1488   long k;
1489
1490   if (x > 0.5 && x < 1.5) { 
1491      xcopy(t1, x - 1);
1492      k = Lg2(t1);
1493   }
1494   else {
1495      k = NumBits(Lg2(x)); 
1496   }
1497
1498   k += Lg2(y);
1499
1500   if (k > NTL_BITS_PER_LONG+10) Error("RR: overflow");
1501
1502   if (k < 0) k = 0;
1503
1504   
1505   RR::SetPrecision(p + k + 10);
1506
1507   t1 = y*log(x);
1508
1509   RR::SetPrecision(p);
1510
1511   t2 = exp(t1);
1512
1513   res = t2;
1514}
1515
1516
1517void ReallyComputePi(RR& res)
1518{
1519   long p = RR::precision();
1520   RR::SetPrecision(p + NumBits(p) + 10);
1521
1522
1523   RR sum1;
1524
1525   RR s, s1, t, t1;
1526
1527   s = 0;
1528   t = 0.5;
1529   t1 = 0.5;
1530
1531   long i;
1532
1533   for (i = 3; ; i+=2) {
1534      add(s1, s, t);
1535      if (s == s1) break;
1536      xcopy(s, s1);
1537      mul(t1, t1, -0.25);
1538      div(t, t1, i);
1539   }
1540
1541   xcopy(sum1, s);
1542
1543
1544   RR g;
1545
1546   inv(g, to_RR(3)); // g = 1/3
1547
1548   s = 0;
1549   xcopy(t, g);
1550   xcopy(t1, g);
1551
1552   sqr(g, g);
1553   negate(g, g); // g = -1/9
1554
1555   for (i = 3; ; i+=2) {
1556      add(s1, s, t);
1557      if (s == s1) break;
1558      xcopy(s, s1);
1559      mul(t1, t1, g);
1560      div(t, t1, i);
1561   }
1562
1563   add(s, s, sum1);
1564   mul(s, s, 4);
1565
1566   RR::SetPrecision(p);
1567   xcopy(res, s);
1568}
1569
1570void ComputePi(RR& res)
1571{
1572   static long prec = 0;
1573   static RR pi;
1574
1575   long p = RR::precision();
1576
1577   if (prec <= p + 10) {
1578      prec = p + 20;
1579      RR::SetPrecision(prec);
1580      ReallyComputePi(pi);
1581      RR::SetPrecision(p);
1582   }
1583
1584   xcopy(res, pi);
1585}
1586
1587
1588
1589void sin(RR& res, const RR& x)
1590{
1591   if (x == 0) {
1592      res = 0;
1593      return;
1594   }
1595
1596   if (Lg2(x) > 1000) 
1597      Error("sin: sorry...argument too large in absolute value");
1598
1599   long p = RR::precision();
1600
1601   RR pi, t1, f;
1602   RR n;
1603
1604
1605   // we want to make x^2 < 3, so that the series for sin(x)
1606   // converges nicely, without any nasty cancellations in the
1607   // first terms of the series.
1608
1609   RR::SetPrecision(p + NumBits(p) + 10);
1610
1611   if (x*x < 3) {
1612      xcopy(f, x);
1613   }
1614   else {
1615
1616      // we want to write x/pi = n + f, |f| < 1/2....
1617      // but we have to do *this* very carefully, so that f is computed
1618      // to precision > p.  I know, this is sick!
1619   
1620      long p1;
1621   
1622      p1 = p + Lg2(x) + 20;
1623   
1624   
1625      for (;;) {
1626         RR::SetPrecision(p1);
1627         ComputePi(pi);
1628         xcopy(t1, x/pi);
1629         xcopy(n, floor(t1));
1630         xcopy(f, t1 - n);
1631         if (f > 0.5) {
1632            n++;
1633            xcopy(f, t1 - n);
1634         }
1635
1636         if (f == 0 || p1 < p - Lg2(f) + Lg2(n) + 10) {
1637            // we don't have enough bits of f...increase p1 and continue
1638
1639            p1 = p1 + max(20, p1/10);
1640         }
1641         else
1642            break;
1643      }
1644
1645      RR::SetPrecision(p + NumBits(p) + 10);
1646      ComputePi(pi);
1647
1648      xcopy(f, pi * f);
1649     
1650      if (n != 0 && n.exponent() == 0) {
1651         // n is odd, so we negate f, which negates sin(f)
1652
1653         xcopy(f, -f);
1654      }
1655
1656   }
1657
1658   // Boy, that was painful, but now its over, and we can simply apply
1659   // the series for sin(f)
1660
1661   RR t2, s, s1, t;
1662   long i;
1663
1664   s = 0;
1665   xcopy(t, f);
1666
1667   for (i = 3; ; i=i+2) {
1668      add(s1, s, t);
1669      if (s == s1) break;
1670      xcopy(s, s1);
1671      mul(t, t, f);
1672      mul(t, t, f);
1673      div(t, t, i-1);
1674      div(t, t, i);
1675      negate(t, t);
1676   }
1677
1678   RR::SetPrecision(p);
1679   
1680   xcopy(res, s);
1681
1682}
1683
1684void cos(RR& res, const RR& x)
1685{
1686   if (x == 0) {
1687      res = 1;
1688      return;
1689   }
1690
1691   if (Lg2(x) > 1000) 
1692      Error("cos: sorry...argument too large in absolute value");
1693
1694   long p = RR::precision();
1695
1696   RR pi, t1, f;
1697   RR n;
1698
1699   // we want to write x/pi = (n+1/2) + f, |f| < 1/2....
1700   // but we have to do *this* very carefully, so that f is computed
1701   // to precision > p.  I know, this is sick!
1702
1703   long p1;
1704
1705   p1 = p + Lg2(x) + 20;
1706
1707
1708   for (;;) {
1709      RR::SetPrecision(p1);
1710      ComputePi(pi);
1711      xcopy(t1, x/pi);
1712      xcopy(n, floor(t1));
1713      xcopy(f, t1 - (n + 0.5));
1714
1715      if (f == 0 || p1 < p - Lg2(f) + Lg2(n) + 10) {
1716         // we don't have enough bits of f...increase p1 and continue
1717
1718         p1 = p1 + max(20, p1/10);
1719      }
1720      else
1721         break;
1722   }
1723
1724   RR::SetPrecision(p + NumBits(p) + 10);
1725   ComputePi(pi);
1726
1727   xcopy(f, pi * f);
1728   
1729   if (n == 0 || n.exponent() != 0) {
1730      // n is even, so we negate f, which negates sin(f)
1731
1732      xcopy(f, -f);
1733   }
1734
1735   // Boy, that was painful, but now its over, and we can simply apply
1736   // the series for sin(f)
1737
1738   RR t2, s, s1, t;
1739   long i;
1740
1741   s = 0;
1742   xcopy(t, f);
1743
1744   for (i = 3; ; i=i+2) {
1745      add(s1, s, t);
1746      if (s == s1) break;
1747      xcopy(s, s1);
1748      mul(t, t, f);
1749      mul(t, t, f);
1750      div(t, t, i-1);
1751      div(t, t, i);
1752      negate(t, t);
1753   }
1754
1755   RR::SetPrecision(p);
1756   
1757   xcopy(res, s);
1758
1759}
1760
1761
1762NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.