7 #include "cl_integer.h"
17 double cl_double_approx (const cl_I& x)
19 // Method: same as cl_I_to_DF().
20 if (eq(x,0)) { return 0.0; }
21 var cl_signean sign = -(cl_signean)minusp(x); // Vorzeichen
22 var cl_I abs_x = (sign==0 ? x : -x);
23 var uintL exp = integer_length(abs_x); // (integer-length x)
24 // NDS zu |x|>0 bilden:
25 var const uintD* MSDptr;
27 I_to_NDS_nocopy(abs_x, MSDptr=,len=,,cl_false,);
28 // MSDptr/len/LSDptr ist die NDS zu x, len>0.
29 // Führende Digits holen: Brauche DF_mant_len+1 Bits, dazu intDsize
30 // Bits (die NDS kann mit bis zu intDsize Nullbits anfangen).
31 // Dann werden diese Bits um (exp mod intDsize) nach rechts geschoben.
32 var uintD msd = msprefnext(MSDptr); // erstes Digit
33 #if (cl_word_size==64)
34 var uint64 msdd = 0; // weitere min(len-1,64/intDsize) Digits
35 #define NEXT_DIGIT(i) \
36 { if (--len == 0) goto ok; \
37 msdd |= (uint64)msprefnext(MSDptr) << (64-(i+1)*intDsize); \
39 DOCONSTTIMES(64/intDsize,NEXT_DIGIT);
42 var uint32 msdd = 0; // weitere min(len-1,32/intDsize) Digits
43 var uint32 msddf = 0; // weitere maximal 32/intDsize Digits
44 #define NEXT_DIGIT(i) \
45 { if (--len == 0) goto ok; \
46 msdd |= (uint32)msprefnext(MSDptr) << (32-(i+1)*intDsize); \
48 DOCONSTTIMES(32/intDsize,NEXT_DIGIT);
50 #define NEXT_DIGIT(i) \
51 { if (--len == 0) goto ok; \
52 msddf |= (uint32)msprefnext(MSDptr) << (32-(i+1)*intDsize); \
54 DOCONSTTIMES(32/intDsize,NEXT_DIGIT);
58 #if (cl_word_size==64)
59 // Die NDS besteht aus msd, msdd und len weiteren Digits.
60 // Das höchste in 2^64*msd+msdd gesetzte Bit ist Bit Nummer
61 // 63 + (exp mod intDsize).
62 var uintL shiftcount = exp % intDsize;
63 var uint64 mant = // führende 64 Bits
66 : (((uint64)msd << (64-shiftcount)) | (msdd >> shiftcount))
68 // Das höchste in mant gesetzte Bit ist Bit Nummer 63.
69 if ( ((mant & bit(62-DF_mant_len)) ==0) // Bit 10 =0 -> abrunden
70 || ( ((mant & (bit(62-DF_mant_len)-1)) ==0) // Bit 10 =1 und Bits 9..0 =0
71 && ((msdd & (bit(shiftcount)-1)) ==0) // und weitere Bits aus msddf =0
72 && (!test_loop_msp(MSDptr,len)) // und alle weiteren Digits =0
73 // round-to-even, je nach Bit 11 :
74 && ((mant & bit(63-DF_mant_len)) ==0)
77 { mant = mant >> (63-DF_mant_len); }
80 { mant = mant >> (63-DF_mant_len);
82 if (mant >= bit(DF_mant_len+1)) // rounding overflow?
83 { mant = mant>>1; exp = exp+1; }
85 union { dfloat eksplicit; double machine_double; } u;
86 if ((sintL)exp > (sintL)(DF_exp_high-DF_exp_mid))
88 ((sint64)sign & bit(63))
89 | ((uint64)(bit(DF_exp_len)-1) << DF_mant_len); // Infinity
93 ((sint64)sign & bit(63)) /* Vorzeichen */
94 | ((uint64)((sintL)exp+DF_exp_mid) << DF_mant_len) /* Exponent */
95 | ((uint64)mant & (bit(DF_mant_len)-1)); /* Mantisse */
97 return u.machine_double;
99 // Die NDS besteht aus msd, msdd, msddf und len weiteren Digits.
100 // Das höchste in 2^64*msd+2^32*msdd+msddf gesetzte Bit ist Bit Nummer
101 // 63 + (exp mod intDsize).
102 var uintL shiftcount = exp % intDsize;
103 var uint32 manthi; // führende 32 Bits
104 var uint32 mantlo; // nächste 32 Bits
106 { manthi = msdd; mantlo = msddf; }
108 { manthi = ((uint32)msd << (32-shiftcount)) | (msdd >> shiftcount);
109 mantlo = (msdd << (32-shiftcount)) | (msddf >> shiftcount);
111 // Das höchste in mant gesetzte Bit ist Bit Nummer 63.
112 if ( ((mantlo & bit(62-DF_mant_len)) ==0) // Bit 10 =0 -> abrunden
113 || ( ((mantlo & (bit(62-DF_mant_len)-1)) ==0) // Bit 10 =1 und Bits 9..0 =0
114 && ((msddf & (bit(shiftcount)-1)) ==0) // und weitere Bits aus msddf =0
115 && (!test_loop_msp(MSDptr,len)) // und alle weiteren Digits =0
116 // round-to-even, je nach Bit 11 :
117 && ((mantlo & bit(63-DF_mant_len)) ==0)
120 { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
121 manthi = manthi >> (63-DF_mant_len);
125 { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
126 manthi = manthi >> (63-DF_mant_len);
130 if (manthi >= bit(DF_mant_len-32+1)) // rounding overflow?
131 { manthi = manthi>>1; exp = exp+1; }
133 union { dfloat eksplicit; double machine_double; } u;
134 if ((sintL)exp > (sintL)(DF_exp_high-DF_exp_mid))
135 { u.eksplicit.semhi =
136 ((sint32)sign & bit(31))
137 | ((uint32)(bit(DF_exp_len)-1) << (DF_mant_len-32)); // Infinity
141 { u.eksplicit.semhi =
142 ((sint32)sign & bit(31)) /* Vorzeichen */
143 | ((uint32)((sintL)exp+DF_exp_mid) << (DF_mant_len-32)) /* Exponent */
144 | ((uint32)manthi & (bit(DF_mant_len-32)-1)); /* Mantisse */
145 u.eksplicit.mlo = mantlo;
147 return u.machine_double;