7 #include "cln/rational.h"
14 #include "cln/integer.h"
20 double double_approx (const cl_RA& x)
22 // Method: same as cl_RA_to_DF().
25 return double_approx(x);
29 union { dfloat eksplicit; double machine_double; } u;
30 var cl_I a = numerator(x); // +/- a
31 var const cl_I& b = denominator(x); // b
32 var cl_signean sign = -(cl_signean)minusp(a); // Vorzeichen
33 if (!(sign==0)) { a = -a; } // Betrag nehmen, liefert a
34 var sintC lendiff = (sintC)integer_length(a) // (integer-length a)
35 - (sintC)integer_length(b); // (integer-length b)
36 if (lendiff > DF_exp_high-DF_exp_mid) // Exponent >= n-m > Obergrenze ?
38 #if (cl_word_size==64)
40 ((sint64)sign & bit(63))
41 | ((uint64)(bit(DF_exp_len)-1) << DF_mant_len); // Infinity
44 ((sint32)sign & bit(31))
45 | ((uint32)(bit(DF_exp_len)-1) << (DF_mant_len-32)); // Infinity
48 return u.machine_double;
50 if (lendiff < DF_exp_low-DF_exp_mid-2) // Exponent <= n-m+2 < Untergrenze ?
52 #if (cl_word_size==64)
53 u.eksplicit = ((sint64)sign & bit(63)); // 0.0
55 u.eksplicit.semhi = ((sint32)sign & bit(31)); // 0.0
58 return u.machine_double;
62 if (lendiff >= DF_mant_len+2)
64 { nenner = ash(b,lendiff - (DF_mant_len+2)); // (ash b n-m-54)
68 { zaehler = ash(a,(DF_mant_len+2) - lendiff); // (ash a -n+m+54)
71 // Division zaehler/nenner durchführen:
72 var cl_I_div_t q_r = cl_divide(zaehler,nenner);
73 var cl_I& q = q_r.quotient;
74 var cl_I& r = q_r.remainder;
75 // 2^53 <= q < 2^55, also ist q Bignum mit ceiling(55/intDsize) Digits.
76 var const uintD* ptr = BN_MSDptr(q);
77 #if (cl_word_size==64)
78 var uint64 mant = get_max64_Dptr(55,ptr);
79 if (mant >= bit(DF_mant_len+2))
80 // 2^54 <= q < 2^55, schiebe um 2 Bits nach rechts
81 { var uint64 rounding_bits = mant & (bit(2)-1);
82 lendiff = lendiff+1; // Exponent := n-m+1
84 if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
85 || ( (rounding_bits == bit(1)) // 10
86 && (eq(r,0)) // und genau halbzahlig (r=0)
87 && ((mant & bit(0)) ==0) // -> round-to-even
96 { var uint64 rounding_bit = mant & bit(0);
98 if ( (rounding_bit == 0) // 0 wird abgerundet
99 || ( (eq(r,0)) // genau halbzahlig (r=0)
100 && ((mant & bit(0)) ==0) // -> round-to-even
110 if (mant >= bit(DF_mant_len+1)) // rounding overflow?
111 { mant = mant>>1; lendiff = lendiff+1; }
114 if (lendiff < (sintL)(DF_exp_low-DF_exp_mid))
115 { u.eksplicit = ((sint64)sign & bit(63)); } // 0.0
116 else if (lendiff > (sintL)(DF_exp_high-DF_exp_mid))
118 ((sint64)sign & bit(63))
119 | ((uint64)(bit(DF_exp_len)-1) << DF_mant_len); // Infinity
123 ((sint64)sign & bit(63)) /* Vorzeichen */
124 | ((uint64)(lendiff+DF_exp_mid) << DF_mant_len) /* Exponent */
125 | ((uint64)mant & (bit(DF_mant_len)-1)); /* Mantisse */
127 return u.machine_double;
129 var uint32 manthi = get_max32_Dptr(23,ptr);
130 var uint32 mantlo = get_32_Dptr(ptr mspop ceiling(23,intDsize));
131 if (manthi >= bit(DF_mant_len-32+2))
132 // 2^54 <= q < 2^55, schiebe um 2 Bits nach rechts
133 { var uintL rounding_bits = mantlo & (bit(2)-1);
134 lendiff = lendiff+1; // Exponent := n-m+1
135 mantlo = (mantlo >> 2) | (manthi << 30); manthi = manthi >> 2;
136 if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
137 || ( (rounding_bits == bit(1)) // 10
138 && (eq(r,0)) // und genau halbzahlig (r=0)
139 && ((mantlo & bit(0)) ==0) // -> round-to-even
148 { var uintL rounding_bit = mantlo & bit(0);
149 mantlo = (mantlo >> 1) | (manthi << 31); manthi = manthi >> 1;
150 if ( (rounding_bit == 0) // 0 wird abgerundet
151 || ( (eq(r,0)) // genau halbzahlig (r=0)
152 && ((mantlo & bit(0)) ==0) // -> round-to-even
164 if (manthi >= bit(DF_mant_len-32+1)) // rounding overflow?
165 { manthi = manthi>>1; lendiff = lendiff+1; }
169 if (lendiff < (sintL)(DF_exp_low-DF_exp_mid))
170 { u.eksplicit.semhi = ((sint32)sign & bit(31)); // 0.0
173 else if (lendiff > (sintL)(DF_exp_high-DF_exp_mid))
174 { u.eksplicit.semhi =
175 ((sint32)sign & bit(31))
176 | ((uint32)(bit(DF_exp_len)-1) << (DF_mant_len-32)); // Infinity
180 { u.eksplicit.semhi =
181 ((sint32)sign & bit(31)) /* Vorzeichen */
182 | ((uint32)(lendiff+DF_exp_mid) << (DF_mant_len-32)) /* Exponent */
183 | ((uint32)manthi & (bit(DF_mant_len-32)-1)); /* Mantisse */
184 u.eksplicit.mlo = mantlo;
186 return u.machine_double;