13 #include "cln/float.h"
14 #include "cln/rational.h"
15 #include "cln/integer.h"
21 // Methode (rekursiv dargestellt):
22 // Falls x rational ist: x.
24 // Falls x<0.0: (- (rationalize (- x)))
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.
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.)
50 // if c>=b then k:=c-1, "Ziffer k", (a,b) := (1/(b-k),1/(a-k)), goto Schleife
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.
57 inline const cl_RA rationalize (const cl_RA& x)
59 // x rational -> x als Ergebnis.
63 inline const cl_RA rationalize (const cl_F& x)
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;
71 var cl_I y = ash(m,e);
72 if (minusp(s)) { y = -y; }
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]
91 var cl_I k = minus1(c); // k = c-1
93 var cl_I p_iplus1 = k * p_i + p_iminus1;
94 p_iminus1 = p_i; p_i = p_iplus1;
97 var cl_I q_iplus1 = k * q_i + q_iminus1;
98 q_iminus1 = q_i; q_i = q_iplus1;
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;
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;
111 return I_I_to_RA(p_last,q_last); // +-p[i] / q[i] bilden
114 const cl_RA rationalize (const cl_R& x)
115 GEN_R_OP1_2(x, rationalize, return)