source: git/ntl/src/RR.c @ f48c17

spielwiese
Last change on this file since f48c17 was f48c17, checked in by Hans Schönemann <hannes@…>, 17 years ago
*hannes: NTL 5.4.1 git-svn-id: file:///usr/local/Singular/svn/trunk@10338 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 27.4 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
397
398void SqrRoot(RR& z, const RR& a)
399{
400   if (sign(a) < 0)
401      Error("RR: attempt to take square root of negative number");
402
403   if (IsZero(a)) {
404      clear(z);
405      return;
406   }
407
408   RR t;
409   ZZ T1, T2;
410   long k;
411
412   k = 2*RR::prec - NumBits(a.x) + 1;
413
414   if (k < 0) k = 0;
415
416   if ((a.e - k) & 1) k++;
417
418   LeftShift(T1, a.x, k);
419   // since k >= 2*prec - bits(a) + 1, T1 has at least 2*prec+1 bits,           
420   // thus T1 >= 2^(2*prec)                                                     
421
422   SqrRoot(t.x, T1); // t.x >= 2^prec thus t.x contains the round bit           
423   t.e = (a.e - k)/2;
424   sqr(T2, t.x); 
425
426   // T1-T2 is the (lower part of the) sticky bit                               
427   normalize(z, t, T2 < T1);
428}
429
430
431
432
433void SqrRootPrec(RR& x, const RR& a, long p)
434{
435   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
436      Error("SqrRootPrec: bad precsion");
437
438   long old_p = RR::prec;
439   RR::prec = p;
440   SqrRoot(x, a);
441   RR::prec = old_p;
442}
443
444
445
446
447void swap(RR& a, RR& b)
448{
449   swap(a.x, b.x);
450   swap(a.e, b.e);
451}
452
453long compare(const RR& a, const RR& b)
454{
455   static RR t;
456
457   SubPrec(t, a, b, 1);
458   return sign(t);
459}
460
461
462
463long operator==(const RR& a, const RR& b) 
464{
465   return a.e == b.e && a.x == b.x;
466}
467
468
469void trunc(RR& z, const RR& a)
470{
471   static RR t;
472
473   if (a.e >= 0) 
474      xcopy(z, a);
475   else {
476      RightShift(t.x, a.x, -a.e);
477      t.e = 0;
478      xcopy(z, t);
479   }
480}
481
482void TruncPrec(RR& x, const RR& a, const RR& b, long p)
483{
484   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
485      Error("TruncPrec: bad precsion");
486
487   long old_p = RR::prec;
488   RR::prec = p;
489   trunc(x, a);
490   RR::prec = old_p;
491}
492
493void floor(RR& z, const RR& a)
494{
495   static RR t;
496
497   if (a.e >= 0) 
498      xcopy(z, a);
499   else {
500      RightShift(t.x, a.x, -a.e);
501      if (sign(a.x) < 0)
502         add(t.x, t.x, -1);
503      t.e = 0;
504      xcopy(z, t);
505   }
506}
507
508void FloorPrec(RR& x, const RR& a, const RR& b, long p)
509{
510   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
511      Error("FloorPrec: bad precsion");
512
513   long old_p = RR::prec;
514   RR::prec = p;
515   floor(x, a);
516   RR::prec = old_p;
517}
518
519void ceil(RR& z, const RR& a)
520{
521   static RR t;
522
523   if (a.e >= 0)
524      xcopy(z, a);
525   else {
526      RightShift(t.x, a.x, -a.e);
527      if (sign(a.x) > 0)
528         add(t.x, t.x, 1);
529      t.e = 0;
530      xcopy(z, t);
531   }
532}
533
534void CeilPrec(RR& x, const RR& a, const RR& b, long p)
535{
536   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
537      Error("CeilPrec: bad precsion");
538
539   long old_p = RR::prec;
540   RR::prec = p;
541   ceil(x, a);
542   RR::prec = old_p;
543}
544
545void round(RR& z, const RR& a)
546{
547   if (a.e >= 0) {
548      xcopy(z, a);
549      return;
550   }
551
552   long len = NumBits(a.x);
553
554   if (-a.e > len) {
555      z = 0;
556      return;
557   }
558
559   if (-a.e == len) {
560      if (len == 1)
561         z = 0;
562      else
563         z = sign(a.x);
564
565      return;
566   }
567
568   static RR t;
569   ConvPrec(t, a, len+a.e);
570   xcopy(z, t);
571}
572
573void RoundPrec(RR& x, const RR& a, const RR& b, long p)
574{
575   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
576      Error("RoundPrec: bad precsion");
577
578   long old_p = RR::prec;
579   RR::prec = p;
580   round(x, a);
581   RR::prec = old_p;
582}
583
584
585   
586
587void conv(RR& z, const ZZ& a)
588{
589   normalize1(z, a, 0, RR::prec, 0);
590}
591
592void ConvPrec(RR& x, const ZZ& a, long p)
593{
594   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
595      Error("ConvPrec: bad precsion");
596
597   long old_p = RR::prec;
598   RR::prec = p;
599   conv(x, a);
600   RR::prec = old_p;
601}
602
603
604void conv(RR& z, long a)
605{
606   if (a == 0) {
607      clear(z);
608      return;
609   }
610
611   if (a == 1) {
612      set(z);
613      return;
614   }
615
616   static ZZ t;
617   t = a;
618   conv(z, t);
619}
620
621void ConvPrec(RR& x, long a, long p)
622{
623   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
624      Error("ConvPrec: bad precsion");
625
626   long old_p = RR::prec;
627   RR::prec = p;
628   conv(x, a);
629   RR::prec = old_p;
630}
631
632void conv(RR& z, unsigned long a)
633{
634   if (a == 0) {
635      clear(z);
636      return;
637   }
638
639   if (a == 1) {
640      set(z);
641      return;
642   }
643
644   static ZZ t;
645   conv(t, a);
646   conv(z, t);
647}
648
649void ConvPrec(RR& x, unsigned long a, long p)
650{
651   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
652      Error("ConvPrec: bad precsion");
653
654   long old_p = RR::prec;
655   RR::prec = p;
656   conv(x, a);
657   RR::prec = old_p;
658}
659
660
661void conv(RR& z, double a)
662{
663   if (a == 0) {
664      clear(z);
665      return;
666   }
667
668   if (a == 1) {
669      set(z);
670      return;
671   }
672
673   if (!IsFinite(&a))
674      Error("RR: conversion of a non-finite double");
675
676   int e;
677   double f;
678   static RR t;
679
680   f = frexp(a, &e);
681
682   f = f * NTL_FDOUBLE_PRECISION;
683   f = f * 4;
684
685   conv(t.x, f);
686   t.e = e - (NTL_DOUBLE_PRECISION + 1);
687
688   xcopy(z, t);
689}
690
691void ConvPrec(RR& x, double a, long p)
692{
693   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
694      Error("ConvPrec: bad precsion");
695
696   long old_p = RR::prec;
697   RR::prec = p;
698   conv(x, a);
699   RR::prec = old_p;
700}
701
702
703void conv(ZZ& z, const RR& a)
704{
705   if (a.e >= 0) 
706      LeftShift(z, a.x, a.e);
707   else {
708      long sgn = sign(a.x);
709      RightShift(z, a.x, -a.e);
710      if (sgn < 0)
711         sub(z, z, 1);
712   }
713}
714
715void CeilToZZ(ZZ& z, const RR& a)
716{
717   if (a.e >= 0)
718      LeftShift(z, a.x, a.e);
719   else {
720      long sgn = sign(a.x);
721      RightShift(z, a.x, -a.e);
722      if (sgn > 0)
723         add(z, z, 1);
724   }
725}
726
727void TruncToZZ(ZZ& z, const RR& a)
728{
729   if (a.e >= 0)
730      LeftShift(z, a.x, a.e);
731   else 
732      RightShift(z, a.x, -a.e);
733}
734
735
736void RoundToZZ(ZZ& z, const RR& a)
737{
738   if (a.e >= 0) {
739      LeftShift(z, a.x, a.e);
740      return;
741   }
742
743   long len = NumBits(a.x);
744
745   if (-a.e > len) {
746      z = 0;
747      return;
748   }
749
750   if (-a.e == len) {
751      if (len == 1)
752         z = 0;
753      else
754         z = sign(a.x);
755
756      return;
757   }
758
759   static RR t;
760
761   ConvPrec(t, a, len+a.e);
762
763   LeftShift(z, t.x, t.e);
764}
765
766
767void conv(long& z, const RR& a)
768{
769   ZZ t;
770   conv(t, a);
771   conv(z, t);
772}
773
774void conv(double& z, const RR& aa)
775{
776   double x;
777   static RR a;
778
779   ConvPrec(a, aa, NTL_DOUBLE_PRECISION);
780   // round to NTL_DOUBLE_PRECISION bits to avoid double overflow
781
782   conv(x, a.x);
783   z = _ntl_ldexp(x, a.e);
784}
785
786
787void add(RR& z, const RR& a, double b)
788{
789   static RR B;
790   B = b;
791   add(z, a, B);
792}
793
794
795
796void sub(RR& z, const RR& a, double b)
797{
798   static RR B;
799   B = b;
800   sub(z, a, B);
801}
802
803void sub(RR& z, double a, const RR& b)
804{
805   static RR A;
806   A = a;
807   sub(z, A, b);
808}
809
810
811
812void mul(RR& z, const RR& a, double b)
813{
814   static RR B;
815   B = b;
816   mul(z, a, B);
817}
818
819
820void div(RR& z, const RR& a, double b)
821{
822   static RR B;
823   B = b;
824   div(z, a, B);
825}
826
827void div(RR& z, double a, const RR& b)
828{
829   static RR A;
830   A = a;
831   div(z, A, b);
832}
833
834
835void inv(RR& z, const RR& a)
836{
837   static RR one = to_RR(1);
838   div(z, one, a);
839}
840
841void InvPrec(RR& x, const RR& a, long p)
842{
843   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
844      Error("InvPrec: bad precsion");
845
846   long old_p = RR::prec;
847   RR::prec = p;
848   inv(x, a);
849   RR::prec = old_p;
850}
851
852
853long compare(const RR& a, double b)
854{
855   if (b == 0) return sign(a);
856
857   static RR B;
858   B = b;
859   return compare(a, B);
860}
861
862
863long operator==(const RR& a, double b) 
864{
865   if (b == 0) return IsZero(a);
866   if (b == 1) return IsOne(a);
867
868   static RR B;
869   B = b;
870   return a == B;
871}
872
873
874void power(RR& z, const RR& a, long e)
875{
876   RR b, res;
877
878   long n = NumBits(e);
879
880   long p = RR::precision();
881   RR::SetPrecision(p + n + 10);
882
883   xcopy(b, a);
884
885   set(res);
886   long i;
887
888   for (i = n-1; i >= 0; i--) {
889      sqr(res, res);
890      if (bit(e, i))
891         mul(res, res, b);
892   }
893
894   RR::SetPrecision(p);
895
896   if (e < 0) 
897      inv(z, res);
898   else
899      xcopy(z, res);
900}
901
902
903void conv(RR& z, const xdouble& a)
904{
905   conv(z, a.mantissa());
906
907   if (a.exponent() >  ((2*NTL_OVFBND)/(2*NTL_XD_HBOUND_LOG))) 
908      Error("RR: overlow");
909
910   if (a.exponent() < -((2*NTL_OVFBND)/(2*NTL_XD_HBOUND_LOG))) 
911      Error("RR: underflow");
912
913   z.e += a.exponent()*(2*NTL_XD_HBOUND_LOG);
914
915   if (z.e >= NTL_OVFBND)
916      Error("RR: overflow");
917
918   if (z.e <= -NTL_OVFBND)
919      Error("RR: underflow");
920}
921
922void ConvPrec(RR& x, const xdouble& a, long p)
923{
924   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
925      Error("ConvPrec: bad precsion");
926
927   long old_p = RR::prec;
928   RR::prec = p;
929   conv(x, a);
930   RR::prec = old_p;
931}
932
933
934
935void conv(xdouble& z, const RR& a)
936{
937   xdouble x;
938   xdouble y;
939
940   conv(x, a.x);
941   power2(y, a.e);
942   z = x*y;
943}
944     
945void power2(RR& z, long e)
946{
947   if (e >= NTL_OVFBND)
948      Error("RR: overflow");
949
950   if (e <= -NTL_OVFBND)
951      Error("RR: underflow");
952
953   set(z.x); 
954   z.e = e;
955}
956
957void conv(RR& z, const quad_float& a)
958{
959   static RR hi, lo, res;
960
961   ConvPrec(hi, a.hi, NTL_DOUBLE_PRECISION);
962   ConvPrec(lo, a.lo, NTL_DOUBLE_PRECISION);
963
964   add(res, hi, lo);
965
966   z = res;
967}
968
969void ConvPrec(RR& x, const quad_float& a, long p)
970{
971   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
972      Error("ConvPrec: bad precsion");
973
974   long old_p = RR::prec;
975   RR::prec = p;
976   conv(x, a);
977   RR::prec = old_p;
978}
979
980
981void conv(quad_float& z, const RR& a)
982{
983   long old_p;
984   static RR a_hi, a_lo;
985
986   old_p = RR::prec;
987
988   ConvPrec(a_hi, a, NTL_DOUBLE_PRECISION);  // high order bits
989   SubPrec(a_lo, a, a_hi, NTL_DOUBLE_PRECISION);  // low order bits
990
991   z = to_quad_float(a_hi.x)*power2_quad_float(a_hi.e) +
992       to_quad_float(a_lo.x)*power2_quad_float(a_lo.e);
993}
994
995void conv(RR& x, const char *s)
996{
997   long c;
998   long cval;
999   long sign;
1000   ZZ a, b;
1001   long i = 0;
1002
1003   if (!s) Error("bad RR input");
1004
1005
1006   c = s[i];
1007   while (IsWhiteSpace(c)) {
1008      i++;
1009      c = s[i];
1010   }
1011
1012   if (c == '-') {
1013      sign = -1;
1014      i++;
1015      c = s[i];
1016   }
1017   else
1018      sign = 1;
1019
1020   long got1 = 0;
1021   long got_dot = 0;
1022   long got2 = 0;
1023
1024   a = 0;
1025   b = 1;
1026
1027   cval = CharToIntVal(c);
1028
1029   if (cval >= 0 && cval <= 9) {
1030      got1 = 1;
1031
1032      while (cval >= 0 && cval <= 9) {
1033         mul(a, a, 10);
1034         add(a, a, cval);
1035         i++;
1036         c = s[i];
1037         cval = CharToIntVal(c);
1038      }
1039   }
1040
1041   if (c == '.') {
1042      got_dot = 1;
1043
1044      i++;
1045      c = s[i];
1046      cval = CharToIntVal(c);
1047
1048      if (cval >= 0 && cval <= 9) {
1049         got2 = 1;
1050   
1051         while (cval >= 0 && cval <= 9) {
1052            mul(a, a, 10);
1053            add(a, a, cval);
1054            mul(b, b, 10);
1055            i++;
1056            c = s[i];
1057            cval = CharToIntVal(c);
1058         }
1059      }
1060   }
1061
1062   if (got_dot && !got1 && !got2)  Error("bad RR input");
1063
1064   ZZ e;
1065
1066   long got_e = 0;
1067   long e_sign;
1068
1069   if (c == 'e' || c == 'E') {
1070      got_e = 1;
1071
1072      i++;
1073      c = s[i];
1074
1075      if (c == '-') {
1076         e_sign = -1;
1077         i++;
1078         c = s[i];
1079      }
1080      else if (c == '+') {
1081         e_sign = 1;
1082         i++;
1083         c = s[i];
1084      }
1085      else
1086         e_sign = 1;
1087
1088
1089      cval = CharToIntVal(c);
1090
1091      if (cval < 0 || cval > 9) Error("bad RR input");
1092
1093      e = 0;
1094      while (cval >= 0 && cval <= 9) {
1095         mul(e, e, 10);
1096         add(e, e, cval);
1097         i++;
1098         c = s[i];
1099         cval = CharToIntVal(c);
1100      }
1101   }
1102
1103   if (!got1 && !got2 && !got_e) Error("bad RR input");
1104
1105   RR t1, t2, v;
1106
1107   long old_p = RR::precision();
1108
1109   if (got1 || got2) {
1110      ConvPrec(t1, a, max(NumBits(a), 1));
1111      ConvPrec(t2, b, NumBits(b));
1112      if (got_e)
1113         RR::SetPrecision(old_p + 10);
1114
1115      div(v, t1, t2);
1116   }
1117   else
1118      set(v);
1119
1120   if (sign < 0)
1121      negate(v, v);
1122
1123   if (got_e) {
1124      if (e >= NTL_OVFBND) Error("RR input overflow");
1125      long E;
1126      conv(E, e);
1127      if (e_sign < 0) E = -E;
1128      RR::SetPrecision(old_p + 10);
1129      power(t1, to_RR(10), E);
1130      mul(v, v, t1);
1131      RR::prec = old_p;
1132   }
1133
1134   xcopy(x, v);
1135}
1136
1137void ConvPrec(RR& x, const char *s, long p)
1138{
1139   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
1140      Error("ConvPrec: bad precsion");
1141
1142   long old_p = RR::prec;
1143   RR::prec = p;
1144   conv(x, s);
1145   RR::prec = old_p;
1146}
1147
1148
1149void ReallyComputeE(RR& res)
1150{
1151   long p = RR::precision();
1152   RR::SetPrecision(p + NumBits(p) + 10);
1153
1154   RR s, s1, t;
1155
1156   s = 1;
1157   t = 1;
1158
1159   long i;
1160
1161   for (i = 2; ; i++) {
1162      add(s1, s, t);
1163      if (s == s1) break;
1164      xcopy(s, s1);
1165      div(t, t, i);
1166   }
1167
1168   RR::SetPrecision(p);
1169   xcopy(res, s);
1170}
1171
1172void ComputeE(RR& res)
1173{
1174   static long prec = 0;
1175   static RR e;
1176
1177   long p = RR::precision();
1178
1179   if (prec <= p + 10) {
1180      prec = p + 20;
1181      RR::SetPrecision(prec);
1182      ReallyComputeE(e);
1183      RR::SetPrecision(p);
1184   }
1185
1186   xcopy(res, e);
1187}
1188
1189
1190void exp(RR& res, const RR& x)
1191{
1192   if (x >= NTL_OVFBND || x <= -NTL_OVFBND)
1193      Error("RR: overflow");
1194
1195   long p = RR::precision();
1196
1197   // step 0: write x = n + f, n an integer and |f| <= 1/2
1198   //    careful -- we want to compute f to > p bits of precision
1199
1200
1201   RR f, nn;
1202   RR::SetPrecision(NTL_BITS_PER_LONG);
1203   round(nn, x);
1204   RR::SetPrecision(p + 10);
1205   sub(f, x, nn);
1206   long n = to_long(nn);
1207
1208   // step 1: calculate t1 = e^n by repeated squaring
1209
1210   RR::SetPrecision(p + NumBits(n) + 10);
1211
1212   RR e;
1213   ComputeE(e);
1214
1215   RR::SetPrecision(p + 10);
1216
1217   RR t1;
1218   power(t1, e, n); 
1219
1220   // step 2: calculate t2 = e^f using Taylor series expansion
1221
1222   RR::SetPrecision(p + NumBits(p) + 10);
1223
1224   RR t2, s, s1, t;
1225   long i;
1226
1227   s = 0;
1228   t = 1;
1229
1230   for (i = 1; ; i++) {
1231      add(s1, s, t);
1232      if (s == s1) break;
1233      xcopy(s, s1);
1234      mul(t, t, f);
1235      div(t, t, i);
1236   }
1237
1238   xcopy(t2, s);
1239
1240   RR::SetPrecision(p);
1241
1242   mul(res, t1, t2);
1243}
1244
1245
1246
1247void ReallyComputeLn2(RR& res)
1248{
1249   long p = RR::precision();
1250   RR::SetPrecision(p + NumBits(p) + 10);
1251
1252   RR s, s1, t, t1;
1253
1254   s = 0;
1255   t = 0.5;
1256   t1 = 0.5;
1257
1258   long i;
1259
1260   for (i = 2; ; i++) {
1261      add(s1, s, t);
1262      if (s == s1) break;
1263      xcopy(s, s1);
1264      mul(t1, t1, 0.5);
1265      div(t, t1, i);
1266   }
1267
1268   RR::SetPrecision(p);
1269   xcopy(res, s);
1270}
1271
1272
1273void ComputeLn2(RR& res)
1274{
1275   static long prec = 0;
1276   static RR ln2;
1277
1278   long p = RR::precision();
1279
1280   if (prec <= p + 10) {
1281      prec = p + 20;
1282      RR::SetPrecision(prec);
1283      ReallyComputeLn2(ln2);
1284      RR::SetPrecision(p);
1285   }
1286
1287   xcopy(res, ln2);
1288}
1289
1290long Lg2(const RR& x)
1291{
1292   return NumBits(x.mantissa()) + x.exponent();
1293}
1294
1295void log(RR& res, const RR& x)
1296{
1297   if (x <= 0) Error("argument to log must be positive");
1298
1299   long p = RR::precision();
1300
1301   RR::SetPrecision(p + NumBits(p) + 10);
1302
1303   RR y;
1304   long n;
1305
1306   // re-write x = 2^n * (1 - y), where -1/2 < y < 1/4  (so 3/4 < 1-y < 3/2)
1307
1308   if (x > 0.75 && x < 1.5) {
1309      n = 0;
1310      sub(y, 1, x);
1311   }
1312   else {
1313      n = Lg2(x) - 1;
1314      RR t;
1315      power2(t, -n);
1316      mul(t, t, x);
1317      while (t > 1.5) {
1318         mul(t, t, 0.5);
1319         n++;
1320      }
1321
1322      sub(y, 1, t);
1323   }
1324
1325   // compute s = - ln(1-y) by power series expansion
1326
1327   RR s, s1, t, t1;
1328
1329   s = 0;
1330   xcopy(t, y);
1331   xcopy(t1, y);
1332
1333   long i;
1334
1335   for (i = 2; ; i++) {
1336      add(s1, s, t);
1337      if (s == s1) break;
1338      xcopy(s, s1);
1339      mul(t1, t1, y);
1340      div(t, t1, i);
1341   }
1342
1343   if (n == 0) 
1344      t = 0;
1345   else {
1346      ComputeLn2(t);
1347      mul(t, t, n);
1348   }
1349
1350   RR::SetPrecision(p);
1351
1352   sub(res, t, s);
1353}
1354
1355
1356void ComputeLn10(RR& res)
1357{
1358   static long prec = 0;
1359   static RR ln10;
1360
1361   long p = RR::precision();
1362
1363   if (prec <= p + 10) {
1364      prec = p + 20;
1365      RR::SetPrecision(prec);
1366      log(ln10, to_RR(10));
1367      RR::SetPrecision(p);
1368   }
1369
1370   xcopy(res, ln10);
1371}
1372
1373void log10(RR& res, const RR& x)
1374{
1375   long p = RR::precision();
1376   RR::SetPrecision(p + 10);
1377
1378   RR ln10, t1, t2;
1379   ComputeLn10(ln10);
1380
1381   log(t1, x);
1382   div(t2, t1, ln10);
1383
1384   RR::SetPrecision(p);
1385
1386   xcopy(res, t2);
1387}
1388
1389
1390void expm1(RR& res, const RR& x)
1391{
1392   long p = RR::precision();
1393
1394   if (x < -0.5 || x > 0.5) {
1395      RR t;
1396      RR::SetPrecision(p + 10);
1397      exp(t, x);
1398      RR::SetPrecision(p);
1399      sub(res, t, 1);
1400      return;
1401   }
1402
1403
1404   RR::SetPrecision(p + NumBits(p) + 10);
1405
1406   RR f;
1407
1408   xcopy(f, x);
1409
1410   RR s, s1, t;
1411   long i;
1412
1413   s = 0;
1414   xcopy(t, f);
1415
1416   for (i = 2; ; i++) {
1417      add(s1, s, t);
1418      if (s == s1) break;
1419      xcopy(s, s1);
1420      mul(t, t, f);
1421      div(t, t, i);
1422   }
1423
1424   RR::SetPrecision(p);
1425
1426   xcopy(res, s);
1427}
1428
1429
1430
1431void log1p(RR& res, const RR& x)
1432{
1433   long p = RR::precision();
1434   RR y;
1435
1436   if (x < -0.5 || x > 0.5) {
1437      RR::SetPrecision(p + 10);
1438      log(y, 1 + x);
1439      RR::SetPrecision(p);
1440      xcopy(res, y);
1441      return;
1442   }
1443
1444
1445   RR::SetPrecision(p + NumBits(p) + 10);
1446
1447
1448   negate(y, x);
1449
1450   // compute s = - ln(1-y) by power series expansion
1451
1452   RR s, s1, t, t1;
1453
1454   s = 0;
1455   xcopy(t, y);
1456   xcopy(t1, y);
1457
1458   long i;
1459
1460   for (i = 2; ; i++) {
1461      add(s1, s, t);
1462      if (s == s1) break;
1463      xcopy(s, s1);
1464      mul(t1, t1, y);
1465      div(t, t1, i);
1466   }
1467
1468   RR::SetPrecision(p);
1469
1470   negate(res, s);
1471
1472}
1473
1474
1475void pow(RR& res, const RR& x, const RR& y)
1476{
1477   if (y == 0) {
1478      res = 1;
1479      return;
1480   }
1481
1482   if (x == 0) {
1483      res = 0;
1484      return;
1485   }
1486
1487   if (x == 1) {
1488      res = 1;
1489      return;
1490   }
1491
1492   if (x < 0) {
1493      Error("pow: sorry...first argument to pow must be nonnegative");
1494   }
1495
1496   long p = RR::precision();
1497
1498   // calculate working precison...one could use p + NTL_BITS_PER_LONG + 10,
1499   // for example, but we want the behaviour to be machine independent.
1500   // so we calculate instead a rough approximation to log |y log(x)|
1501
1502   RR t1, t2;
1503   long k;
1504
1505   if (x > 0.5 && x < 1.5) { 
1506      xcopy(t1, x - 1);
1507      k = Lg2(t1);
1508   }
1509   else {
1510      k = NumBits(Lg2(x)); 
1511   }
1512
1513   k += Lg2(y);
1514
1515   if (k > NTL_BITS_PER_LONG+10) Error("RR: overflow");
1516
1517   if (k < 0) k = 0;
1518
1519   
1520   RR::SetPrecision(p + k + 10);
1521
1522   t1 = y*log(x);
1523
1524   RR::SetPrecision(p);
1525
1526   t2 = exp(t1);
1527
1528   res = t2;
1529}
1530
1531
1532void ReallyComputePi(RR& res)
1533{
1534   long p = RR::precision();
1535   RR::SetPrecision(p + NumBits(p) + 10);
1536
1537
1538   RR sum1;
1539
1540   RR s, s1, t, t1;
1541
1542   s = 0;
1543   t = 0.5;
1544   t1 = 0.5;
1545
1546   long i;
1547
1548   for (i = 3; ; i+=2) {
1549      add(s1, s, t);
1550      if (s == s1) break;
1551      xcopy(s, s1);
1552      mul(t1, t1, -0.25);
1553      div(t, t1, i);
1554   }
1555
1556   xcopy(sum1, s);
1557
1558
1559   RR g;
1560
1561   inv(g, to_RR(3)); // g = 1/3
1562
1563   s = 0;
1564   xcopy(t, g);
1565   xcopy(t1, g);
1566
1567   sqr(g, g);
1568   negate(g, g); // g = -1/9
1569
1570   for (i = 3; ; i+=2) {
1571      add(s1, s, t);
1572      if (s == s1) break;
1573      xcopy(s, s1);
1574      mul(t1, t1, g);
1575      div(t, t1, i);
1576   }
1577
1578   add(s, s, sum1);
1579   mul(s, s, 4);
1580
1581   RR::SetPrecision(p);
1582   xcopy(res, s);
1583}
1584
1585void ComputePi(RR& res)
1586{
1587   static long prec = 0;
1588   static RR pi;
1589
1590   long p = RR::precision();
1591
1592   if (prec <= p + 10) {
1593      prec = p + 20;
1594      RR::SetPrecision(prec);
1595      ReallyComputePi(pi);
1596      RR::SetPrecision(p);
1597   }
1598
1599   xcopy(res, pi);
1600}
1601
1602
1603
1604void sin(RR& res, const RR& x)
1605{
1606   if (x == 0) {
1607      res = 0;
1608      return;
1609   }
1610
1611   if (Lg2(x) > 1000) 
1612      Error("sin: sorry...argument too large in absolute value");
1613
1614   long p = RR::precision();
1615
1616   RR pi, t1, f;
1617   RR n;
1618
1619
1620   // we want to make x^2 < 3, so that the series for sin(x)
1621   // converges nicely, without any nasty cancellations in the
1622   // first terms of the series.
1623
1624   RR::SetPrecision(p + NumBits(p) + 10);
1625
1626   if (x*x < 3) {
1627      xcopy(f, x);
1628   }
1629   else {
1630
1631      // we want to write x/pi = n + f, |f| < 1/2....
1632      // but we have to do *this* very carefully, so that f is computed
1633      // to precision > p.  I know, this is sick!
1634   
1635      long p1;
1636   
1637      p1 = p + Lg2(x) + 20;
1638   
1639   
1640      for (;;) {
1641         RR::SetPrecision(p1);
1642         ComputePi(pi);
1643         xcopy(t1, x/pi);
1644         xcopy(n, floor(t1));
1645         xcopy(f, t1 - n);
1646         if (f > 0.5) {
1647            n++;
1648            xcopy(f, t1 - n);
1649         }
1650
1651         if (f == 0 || p1 < p - Lg2(f) + Lg2(n) + 10) {
1652            // we don't have enough bits of f...increase p1 and continue
1653
1654            p1 = p1 + max(20, p1/10);
1655         }
1656         else
1657            break;
1658      }
1659
1660      RR::SetPrecision(p + NumBits(p) + 10);
1661      ComputePi(pi);
1662
1663      xcopy(f, pi * f);
1664     
1665      if (n != 0 && n.exponent() == 0) {
1666         // n is odd, so we negate f, which negates sin(f)
1667
1668         xcopy(f, -f);
1669      }
1670
1671   }
1672
1673   // Boy, that was painful, but now its over, and we can simply apply
1674   // the series for sin(f)
1675
1676   RR t2, s, s1, t;
1677   long i;
1678
1679   s = 0;
1680   xcopy(t, f);
1681
1682   for (i = 3; ; i=i+2) {
1683      add(s1, s, t);
1684      if (s == s1) break;
1685      xcopy(s, s1);
1686      mul(t, t, f);
1687      mul(t, t, f);
1688      div(t, t, i-1);
1689      div(t, t, i);
1690      negate(t, t);
1691   }
1692
1693   RR::SetPrecision(p);
1694   
1695   xcopy(res, s);
1696
1697}
1698
1699void cos(RR& res, const RR& x)
1700{
1701   if (x == 0) {
1702      res = 1;
1703      return;
1704   }
1705
1706   if (Lg2(x) > 1000) 
1707      Error("cos: sorry...argument too large in absolute value");
1708
1709   long p = RR::precision();
1710
1711   RR pi, t1, f;
1712   RR n;
1713
1714   // we want to write x/pi = (n+1/2) + f, |f| < 1/2....
1715   // but we have to do *this* very carefully, so that f is computed
1716   // to precision > p.  I know, this is sick!
1717
1718   long p1;
1719
1720   p1 = p + Lg2(x) + 20;
1721
1722
1723   for (;;) {
1724      RR::SetPrecision(p1);
1725      ComputePi(pi);
1726      xcopy(t1, x/pi);
1727      xcopy(n, floor(t1));
1728      xcopy(f, t1 - (n + 0.5));
1729
1730      if (f == 0 || p1 < p - Lg2(f) + Lg2(n) + 10) {
1731         // we don't have enough bits of f...increase p1 and continue
1732
1733         p1 = p1 + max(20, p1/10);
1734      }
1735      else
1736         break;
1737   }
1738
1739   RR::SetPrecision(p + NumBits(p) + 10);
1740   ComputePi(pi);
1741
1742   xcopy(f, pi * f);
1743   
1744   if (n == 0 || n.exponent() != 0) {
1745      // n is even, so we negate f, which negates sin(f)
1746
1747      xcopy(f, -f);
1748   }
1749
1750   // Boy, that was painful, but now its over, and we can simply apply
1751   // the series for sin(f)
1752
1753   RR t2, s, s1, t;
1754   long i;
1755
1756   s = 0;
1757   xcopy(t, f);
1758
1759   for (i = 3; ; i=i+2) {
1760      add(s1, s, t);
1761      if (s == s1) break;
1762      xcopy(s, s1);
1763      mul(t, t, f);
1764      mul(t, t, f);
1765      div(t, t, i-1);
1766      div(t, t, i);
1767      negate(t, t);
1768   }
1769
1770   RR::SetPrecision(p);
1771   
1772   xcopy(res, s);
1773
1774}
1775
1776
1777NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.