src/real/misc/cl_R_rationalize.cc

   1 // rationalize().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_real.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_R.h"
  13 #include "cl_float.h"
  14 #include "cl_rational.h"
  15 #include "cl_integer.h"
  16 #include "cl_RA.h"
  17 #include "cl_I.h"
  18
  19 // Methode (rekursiv dargestellt):
  20 // Falls x rational ist: x.
  21 // Falls x=0.0: 0.
  22 // Falls x<0.0: (- (rationalize (- x)))
  23 // Falls x>0.0:
  24 //   (Integer-Decode-Float x) liefert m,e,s=1.
  25 //   Falls e>=0 : Liefere x=m*2^e als Ergebnis.
  26 //   Suche rationale Zahl zwischen a=(m-1/2)*2^e und b=(m+1/2)*2^e mit
  27 //   möglichst kleinem Zähler und Nenner. (a,b einschließlich, aber da a,b
  28 //   den Nenner 2^(|e|+1) haben, während x selbst den Nenner <=2^|e| hat,
  29 //   können weder a noch b als Ergebnis herauskommen.)
  30 //   Suche also bei gegebenem a,b (0<a<b) Bruch y mit a <= y <= b.
  31 //   Rekursiv:
  32 //     c:=(ceiling a)
  33 //     if c<b then return c      ; weil a<=c<b, c ganz
  34 //            else ; a nicht ganz (sonst c=a<b)
  35 //              k:=c-1 ; k=floor(a), k < a < b <= k+1
  36 //              return y = k + 1/(Bruch zwischen 1/(b-k) und 1/(a-k))
  37 //                                ; wobei 1 <= 1/(b-k) < 1/(a-k)
  38 // Man sieht, daß hierbei eine Kettenbruchentwicklung auftritt.
  39 // Methode (iterativ):
  40 // Falls x rational: x.
  41 // (Integer-Decode-Float x) liefert m,e,s.
  42 // e>=0 -> m*2^e*s als Ergebnis (darin ist x=0.0 inbegriffen).
  43 // Bilde a:=(2*m-1)*2^(e-1) und b:=(2*m+1)*2^(e-1), rationale Zahlen >0,
  44 //   (unkürzbar, da Nenner Zweierpotenz und Zähler ungerade).
  45 // Starte Kettenbruchentwicklung (d.h. p[-1]:=0, p[0]:=1, q[-1]:=1, q[0]:=0, i:=0.)
  46 // Schleife:
  47 //   c:=(ceiling a)
  48 //   if c>=b then k:=c-1, "Ziffer k", (a,b) := (1/(b-k),1/(a-k)), goto Schleife
  49 // "Ziffer c".
  50 // (Dabei bedeutet "Ziffer a" die Iteration
  51 //   i:=i+1, p[i]:=a*p[i-1]+p[i-2], q[i]:=a*q[i-1]+q[i-2].)
  52 // Ende, liefere s * (p[i]/q[i]), das ist wegen der Invarianten
  53 //   p[i]*q[i-1]-p[i-1]*q[i]=(-1)^i  ein bereits gekürzter Bruch.
  54
  55 inline const cl_RA rationalize (const cl_RA& x)
  56 {
  57         // x rational -> x als Ergebnis.
  58         return x;
  59 }
  60
  61 inline const cl_RA rationalize (const cl_F& x)
  62 {
  63         var cl_idecoded_float x_decoded = integer_decode_float(x);
  64         var cl_I& m = x_decoded.mantissa;
  65         var cl_I& e = x_decoded.exponent;
  66         var cl_I& s = x_decoded.sign;
  67         if (!minusp(e)) {
  68                 // e>=0.
  69                 var cl_I y = ash(m,e);
  70                 if (minusp(s)) { y = -y; }
  71                 return y;
  72         }
  73         // e<0.
  74         var cl_I m2 = ash(m,1); // 2*m
  75         var cl_I num1 = minus1(m2); // 2*m-1
  76         var cl_I num2 = plus1(m2); // 2*m+1
  77         var cl_I den = ash(1,plus1(-e)); // 2^(1-e)
  78         var cl_RA a = I_I_to_RT(num1,den); // a := (2*m-1)/(2^(1-e))
  79         var cl_RA b = I_I_to_RT(num2,den); // b := (2*m+1)/(2^(1-e))
  80         var cl_I p_iminus1 = 0;         // p[i-1]
  81         var cl_I p_i       = 1;         // p[i]
  82         var cl_I q_iminus1 = 1;         // q[i-1]
  83         var cl_I q_i       = 0;         // q[i]
  84         var cl_I c;
  85         for (;;) {
  86                 c = ceiling1(a);
  87                 if (c < b)
  88                         break;
  89                 var cl_I k = minus1(c); // k = c-1
  90                 {
  91                         var cl_I p_iplus1 = k * p_i + p_iminus1;
  92                         p_iminus1 = p_i; p_i = p_iplus1;
  93                 }
  94                 {
  95                         var cl_I q_iplus1 = k * q_i + q_iminus1;
  96                         q_iminus1 = q_i; q_i = q_iplus1;
  97                 }
  98                 {
  99                         var cl_RA new_b = recip(a-k); // 1/(a-k)
 100                         var cl_RA new_a = recip(b-k); // 1/(b-k)
 101                         a = new_a; b = new_b;
 102                 }
 103         }
 104         // letzte "Ziffer" k=c :
 105         var cl_I p_last = c * p_i + p_iminus1;
 106         var cl_I q_last = c * q_i + q_iminus1;
 107         if (minusp(s))
 108                 p_last = - p_last;
 109         return I_I_to_RA(p_last,q_last); // +-p[i] / q[i] bilden
 110 }
 111
 112 const cl_RA rationalize (const cl_R& x)
 113 GEN_R_OP1_2(x, rationalize, return)