15 const cl_SF sqrt (const cl_SF& x)
18 // x = 0.0 -> Ergebnis 0.0
19 // Ergebnis-Vorzeichen := positiv,
20 // Ergebnis-Exponent := ceiling(e/2),
22 // Bilde aus [1,m15,...,m0,(19 Nullbits)] bei geradem e,
23 // aus [0,1,m15,...,m0,(18 Nullbits)] bei ungeradem e
24 // die Ganzzahl-Wurzel, eine 18-Bit-Zahl mit einer führenden 1.
25 // Runde das letzte Bit weg:
26 // Bit 0 = 0 -> abrunden,
27 // Bit 0 = 1 und Wurzel exakt -> round-to-even,
28 // Bit 0 = 1 und Rest >0 -> aufrunden.
29 // Dabei um ein Bit nach rechts schieben.
30 // Bei Aufrundung auf 2^17 (rounding overflow) Mantisse um 1 Bit nach rechts
31 // schieben und Exponent incrementieren.
35 SF_decode(x, { return x; }, ,exp=,mant=);
36 // Um die 64-Bit-Ganzzahl-Wurzel ausnutzen zu können, fügen wir beim
37 // Radikanden 46 bzw. 47 statt 18 bzw. 19 Nullbits an.
40 { mant = mant << (31-(SF_mant_len+1)); exp = exp+1; }
43 { mant = mant << (32-(SF_mant_len+1)); }
44 exp = exp >> 1; // exp := exp/2
46 isqrt_64_32(mant,0, mant=,exactp=); // mant := isqrt(mant*2^32), eine 32-Bit-Zahl
47 // Die hinteren 31-SF_mant_len Bits wegrunden:
48 if ( ((mant & bit(30-SF_mant_len)) ==0) // Bit 14 =0 -> abrunden
49 || ( ((mant & (bit(30-SF_mant_len)-1)) ==0) // Bit 14 =1 und Bits 13..0 >0 -> aufrunden
50 && exactp // Bit 14 =1 und Bits 13..0 =0, aber Rest -> aufrunden
51 // round-to-even, je nach Bit 15 :
52 && ((mant & bit(31-SF_mant_len)) ==0)
55 { mant = mant >> (31-SF_mant_len); }
58 { mant = mant >> (31-SF_mant_len);
60 if (mant >= bit(SF_mant_len+1)) // rounding overflow?
61 { mant = mant>>1; exp = exp+1; }
63 return encode_SF(0,exp,mant);