12 #include "cl_LF_impl.h"
14 #include "cln/integer.h"
20 const cl_LF cl_RA_to_LF (const cl_RA& x, uintC len)
24 // x = +/- a/b mit Integers a,b>0:
25 // Sei k,m so gewählt, daß
26 // 2^(k-1) <= a < 2^k, 2^(m-1) <= b < 2^m.
27 // Dann ist 2^(k-m-1) < a/b < 2^(k-m+1).
28 // Ergebnis-Vorzeichen := Vorzeichen von x.
29 // Berechne k=(integer-length a) und m=(integer-length b).
30 // Ergebnis-Exponent := k-m.
32 // Berechne floor(2^(-k+m+16n+1)*a/b) :
33 // Bei k-m>=16n+1 dividiere a durch (ash b (k-m-16n-1)),
34 // bei k-m<16n+1 dividiere (ash a (-k+m+16n+1)) durch b.
35 // Der erste Wert ist >=2^16n, <2^(16n+2).
36 // Falls er >=2^(16n+1) ist, erhöhe Exponent um 1,
37 // runde 2 Bits weg und schiebe dabei um 2 Bits nach rechts;
38 // falls er <2^(16n+1) ist,
39 // runde 1 Bit weg und schiebe dabei um 1 Bit nach rechts.
40 // NB: Wenn a und b länger sind als len, ist dieser Algorithmus weniger
41 // effizient, als cl_float(a,len)/cl_float(b,len) zu berechnen. Aber
42 // es ist wichtig, dass cl_RA_to_LF nicht mehr als 0.5 ulp Fehler hat,
43 // deswegen belassen wir es beim ineffizienten aber exakten Algorithmus.
44 // Wenn es auf Rundungsfehler nicht ankommt, muss der Aufrufer im Fall
45 // ceiling(integer_length(a),intDsize) >= len
46 // && ceiling(integer_length(b),intDsize) >= len
47 // einen anderen Algorithmus wählen.
50 return cl_I_to_LF(x,len);
54 var cl_I a = numerator(x); // +/- a
55 var const cl_I& b = denominator(x); // b
56 var cl_signean sign = -(cl_signean)minusp(a); // Vorzeichen
57 if (!(sign==0)) { a = -a; } // Betrag nehmen, liefert a
58 var sintL lendiff = (sintL)integer_length(a) // (integer-length a)
59 - (sintL)integer_length(b); // (integer-length b)
60 // |lendiff| < intDsize*2^intCsize. Da für LF-Exponenten ein sintL zur
61 // Verfügung steht, braucht man keinen Test auf Overflow oder Underflow.
62 var uintL difflimit = intDsize*(uintL)len + 1; // 16n+1
65 if (lendiff > (sintL)difflimit)
66 // 0 <= k-m-16n-1 < k < intDsize*2^intCsize
67 { nenner = ash(b,(uintL)(lendiff - difflimit));
71 // 0 < -k+m+16n+1 <= m+1 + 16n < intDsize*2^intCsize + intDsize*2^intCsize
72 { zaehler = ash(a,(uintL)(difflimit - lendiff)); // (ash a -k+m+16n+1)
75 // Division zaehler/nenner durchführen:
76 var cl_I_div_t q_r = cl_divide(zaehler,nenner);
77 var cl_I& q = q_r.quotient;
78 var cl_I& r = q_r.remainder;
79 // 2^16n <= q < 2^(16n+2), also ist q Bignum mit n+1 Digits.
80 var Lfloat y = allocate_lfloat(len,lendiff+LF_exp_mid,sign); // neues Long-Float
81 var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
82 {var uintD* q_MSDptr = arrayMSDptr(TheBignum(q)->data,len+1);
83 if (mspref(q_MSDptr,0) == 1) // erstes Digit =1 oder =2,3 ?
84 // 2^16n <= q < 2^(16n+1), also 2^(k-m-1) < a/b < 2^(k-m).
85 { // Mantisse mit einer Schiebeschleife um 1 Bit nach rechts füllen:
86 var uintD rounding_bit =
87 shiftrightcopy_loop_msp(q_MSDptr mspop 1,y_mantMSDptr,len,1,1);
88 if ( (rounding_bit == 0) // herausgeschobenes Bit =0 -> abrunden
89 || ( eq(r,0) // =1 und Rest r > 0 -> aufrunden
91 && ((mspref(y_mantMSDptr,len-1) & bit(0)) ==0)
95 goto auf; // aufrunden
98 // 2^(16n+1) <= q < 2^(16n+2), also 2^(k-m) < a/b < 2^(k-m+1).
99 { // Mantisse mit einer Schiebeschleife um 2 Bit nach rechts füllen:
100 var uintD rounding_bits =
101 shiftrightcopy_loop_msp(q_MSDptr mspop 1,y_mantMSDptr,len,2,mspref(q_MSDptr,0));
102 (TheLfloat(y)->expo)++; // Exponenten incrementieren auf k-m+1
103 if ( ((sintD)rounding_bits >= 0) // herausgeschobenes Bit =0 -> abrunden
104 || ( ((rounding_bits & bit(intDsize-2)) ==0) // =1 und nächstes Bit =1 oder Rest r > 0 -> aufrunden
107 && ((mspref(y_mantMSDptr,len-1) & bit(0)) ==0)
111 goto auf; // aufrunden
115 { if ( inc_loop_lsp(y_mantMSDptr mspop len,len) )
116 // Übertrag durchs Aufrunden
117 { mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
118 (TheLfloat(y)->expo)++; // Exponenten incrementieren
124 // Timings on an i486 33 MHz, running Linux, in 0.01 sec.
125 // First timing: cl_I_to_LF(numerator,len)/cl_I_to_LF(denominator,len)
126 // Second timing: cl_RA_to_LF(x,len)
128 // num_length 50 70 100 200 500
131 // 50 1.86 0.97 1.84 0.97 1.85 0.96 1.86 1.86 1.85 7.14
133 // 70 1.86 1.33 1.85 1.31 1.85 1.32 1.84 1.84 1.85 7.13
135 // 100 1.85 1.85 1.86 1.85 1.85 1.84 1.84 1.84 1.86 7.13
137 // 200 1.85 3.61 1.84 3.61 1.85 3.59 1.85 3.59 1.87 7.12
139 // 500 1.84 7.44 1.84 7.55 1.85 7.56 1.84 7.66 1.86 7.63
141 // We see that cl_RA_to_LF is faster only if
142 // num_length < 2*len && den_length < len
143 // whereas cl_I_to_LF(numerator,len)/cl_I_to_LF(denominator,len) is faster if
144 // num_length > 2*len || den_length > len