7 #include "cln/dfloat.h"
19 const cl_DF sqrt (const cl_DF& x)
22 // x = 0.0 -> Ergebnis 0.0
23 // Ergebnis-Vorzeichen := positiv,
24 // Ergebnis-Exponent := ceiling(e/2),
26 // Bilde aus [1,m51,...,m0,(55 Nullbits)] bei geradem e,
27 // aus [0,1,m51,...,m0,(54 Nullbits)] bei ungeradem e
28 // die Ganzzahl-Wurzel, eine 54-Bit-Zahl mit einer führenden 1.
29 // Runde das letzte Bit weg:
30 // Bit 0 = 0 -> abrunden,
31 // Bit 0 = 1 und Wurzel exakt -> round-to-even,
32 // Bit 0 = 1 und Rest >0 -> aufrunden.
33 // Dabei um ein Bit nach rechts schieben.
34 // Bei Aufrundung auf 2^53 (rounding overflow) Mantisse um 1 Bit nach rechts
35 // schieben und Exponent incrementieren.
36 #if (cl_word_size==64)
40 DF_decode(x, { return x; }, ,exp=,mantx=);
41 // Um die 128-Bit-Ganzzahl-Wurzel ausnutzen zu können, fügen wir beim
42 // Radikanden 74 bzw. 75 statt 54 bzw. 55 Nullbits an.
45 { mantx = mantx << (63-(DF_mant_len+1)); exp = exp+1; }
48 { mantx = mantx << (64-(DF_mant_len+1)); }
49 exp = exp >> 1; // exp := exp/2
50 var uintD mant [128/intDsize];
52 arrayLSref(mant,128/intDsize,1) = mantx;
53 arrayLSref(mant,128/intDsize,0) = 0;
54 #else // (intDsize<=32)
55 set_32_Dptr(arrayMSDptr(mant,128/intDsize),(uint32)(mantx>>32));
56 set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 32/intDsize,(uint32)mantx);
57 set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 2*32/intDsize,0);
58 set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 3*32/intDsize,0);
62 var cl_boolean exactp;
63 UDS_sqrt(arrayMSDptr(mant,128/intDsize),128/intDsize,arrayLSDptr(mant,128/intDsize), &wurzel, exactp=);
64 // wurzel = isqrt(2^74_75 * mant), eine 64-Bit-Zahl.
65 mantx = get_64_Dptr(wurzel.MSDptr);
66 // Die hinteren 63-DF_mant_len Bits wegrunden:
67 if ( ((mantx & bit(62-DF_mant_len)) ==0) // Bit 10 =0 -> abrunden
68 || ( ((mantx & (bit(62-DF_mant_len)-1)) ==0) // Bit 10 =1 und Bits 9..0 >0 -> aufrunden
69 && exactp // Bit 10 =1 und Bits 9..0 =0, aber Rest -> aufrunden
70 // round-to-even, je nach Bit 11 :
71 && ((mantx & bit(63-DF_mant_len)) ==0)
74 { mantx = mantx >> (63-DF_mant_len); }
77 { mantx = mantx >> (63-DF_mant_len);
79 if (mantx >= bit(DF_mant_len+1)) // rounding overflow?
80 { mantx = mantx>>1; exp = exp+1; }
83 return encode_DF(0,exp,mantx);
89 DF_decode2(x, { return x; }, ,exp=,manthi=,mantlo=);
90 // Um die 128-Bit-Ganzzahl-Wurzel ausnutzen zu können, fügen wir beim
91 // Radikanden 74 bzw. 75 statt 54 bzw. 55 Nullbits an.
94 { manthi = (manthi << (63-(DF_mant_len+1))) | (mantlo >> ((DF_mant_len+1)-31));
95 mantlo = mantlo << (63-(DF_mant_len+1));
100 { manthi = (manthi << (64-(DF_mant_len+1))) | (mantlo >> ((DF_mant_len+1)-32));
101 mantlo = mantlo << (64-(DF_mant_len+1));
103 exp = exp >> 1; // exp := exp/2
104 var uintD mant [128/intDsize];
105 #if (intDsize==32) || (intDsize==16) || (intDsize==8)
106 set_32_Dptr(arrayMSDptr(mant,128/intDsize),manthi);
107 set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 32/intDsize,mantlo);
108 set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 2*32/intDsize,0);
109 set_32_Dptr(arrayMSDptr(mant,128/intDsize) mspop 3*32/intDsize,0);
112 ptr = arrayLSDptr(mant,128/intDsize);
113 doconsttimes(64/intDsize, { lsprefnext(ptr) = 0; } );
114 doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)mantlo; mantlo = mantlo>>intDsize; } );
115 doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)manthi; manthi = manthi>>intDsize; } );
120 var cl_boolean exactp;
121 UDS_sqrt(arrayMSDptr(mant,128/intDsize),128/intDsize,arrayLSDptr(mant,128/intDsize), &wurzel, exactp=);
122 // wurzel = isqrt(2^74_75 * mant), eine 64-Bit-Zahl.
123 {var uintD* ptr = wurzel.MSDptr;
124 manthi = get_32_Dptr(ptr); mantlo = get_32_Dptr(ptr mspop 32/intDsize);
126 // Die hinteren 63-DF_mant_len Bits wegrunden:
127 if ( ((mantlo & bit(62-DF_mant_len)) ==0) // Bit 10 =0 -> abrunden
128 || ( ((mantlo & (bit(62-DF_mant_len)-1)) ==0) // Bit 10 =1 und Bits 9..0 >0 -> aufrunden
129 && exactp // Bit 10 =1 und Bits 9..0 =0, aber Rest -> aufrunden
130 // round-to-even, je nach Bit 11 :
131 && ((mantlo & bit(63-DF_mant_len)) ==0)
134 { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
135 manthi = manthi >> (63-DF_mant_len);
139 { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
140 manthi = manthi >> (63-DF_mant_len);
144 if (manthi >= bit(DF_mant_len-32+1)) // rounding overflow?
145 { manthi = manthi>>1; exp = exp+1; }
148 return encode_DF(0,exp,manthi,mantlo);