src/real/misc/cl_R_rationalize.cc

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