1 // float_approx().
2
3 // General includes.
4 #include "base/cl_sysdep.h"
5
6 // Specification.
7 #include "cln/rational.h"
8
9
10 // Implementation.
11
12 #include "float/ffloat/cl_FF.h"
13 #include "rational/cl_RA.h"
14 #include "cln/integer.h"
15 #include "integer/cl_I.h"
16 #include "float/cl_F.h"
17
18 namespace cln {
19
float_approx(const cl_RA & x)20 float float_approx (const cl_RA& x)
21 {
22 // Method: same as cl_RA_to_FF().
23 if (integerp(x)) {
24 DeclareType(cl_I,x);
25 return float_approx(x);
26 }
27 { // x Ratio
28 DeclareType(cl_RT,x);
29 union { ffloat eksplicit; float machine_float; } 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 > FF_exp_high-FF_exp_mid) // Exponent >= n-m > Obergrenze ?
37 { u.eksplicit = make_FF_word(sign,bit(FF_exp_len)-1,0); // Infinity
38 return u.machine_float;
39 }
40 if (lendiff < FF_exp_low-FF_exp_mid-2) // Exponent <= n-m+2 < Untergrenze ?
41 { u.eksplicit = make_FF_word(sign,0,0); // 0.0
42 return u.machine_float;
43 }
44 var cl_I zaehler;
45 var cl_I nenner;
46 if (lendiff >= FF_mant_len+2)
47 // n-m-25>=0
48 { nenner = ash(b,lendiff - (FF_mant_len+2)); // (ash b n-m-25)
49 zaehler = a; // a
50 }
51 else
52 { zaehler = ash(a,(FF_mant_len+2) - lendiff); // (ash a -n+m+25)
53 nenner = b; // b
54 }
55 // Division zaehler/nenner durchführen:
56 var cl_I_div_t q_r = cl_divide(zaehler,nenner);
57 var cl_I& q = q_r.quotient;
58 var cl_I& r = q_r.remainder;
59 // 2^24 <= q < 2^26, also ist q Fixnum oder Bignum mit bn_minlength Digits.
60 var uint32 mant = ((FF_mant_len+3 < cl_value_len)
61 ? FN_to_UV(q)
62 : cl_I_to_UL(q)
63 );
64 if (mant >= bit(FF_mant_len+2))
65 // 2^25 <= q < 2^26, schiebe um 2 Bits nach rechts
66 { var uintL rounding_bits = mant & (bit(2)-1);
67 lendiff = lendiff+1; // Exponent := n-m+1
68 mant = mant >> 2;
69 if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
70 || ( (rounding_bits == bit(1)) // 10
71 && (eq(r,0)) // und genau halbzahlig (r=0)
72 && ((mant & bit(0)) ==0) // -> round-to-even
73 ) )
74 // abrunden
75 goto ab;
76 else
77 // aufrunden
78 goto auf;
79 }
80 else
81 { var uintL rounding_bit = mant & bit(0);
82 mant = mant >> 1;
83 if ( (rounding_bit == 0) // 0 wird abgerundet
84 || ( (eq(r,0)) // genau halbzahlig (r=0)
85 && ((mant & bit(0)) ==0) // -> round-to-even
86 ) )
87 // abrunden
88 goto ab;
89 else
90 // aufrunden
91 goto auf;
92 }
93 auf:
94 mant += 1;
95 if (mant >= bit(FF_mant_len+1)) // rounding overflow?
96 { mant = mant>>1; lendiff = lendiff+1; }
97 ab:
98 // Fertig.
99 if (lendiff < (sintL)(FF_exp_low-FF_exp_mid))
100 { u.eksplicit = make_FF_word(sign,0,0); }
101 else if (lendiff > (sintL)(FF_exp_high-FF_exp_mid))
102 { u.eksplicit = make_FF_word(sign,bit(FF_exp_len)-1,0); } // Infinity
103 else
104 { u.eksplicit = make_FF_word(sign,lendiff+FF_exp_mid,mant); }
105 return u.machine_float;
106 }}
107
108 } // namespace cln
109