/cln-1.3.2/src/float/ffloat/algebraic/cl_FF_sqrt.cc

 1// sqrt().
 2
 3// General includes.
 4#include "base/cl_sysdep.h"
 5
 6// Specification.
 7#include "cln/ffloat.h"
 8
 9
10// Implementation.
11
12#include "float/ffloat/cl_FF.h"
13#include "float/cl_F.h"
14#include "base/cl_low.h"
15
16namespace cln {
17
18const cl_FF sqrt (const cl_FF& x)
19{
20// Methode:
21// x = 0.0 -> Ergebnis 0.0
22// Ergebnis-Vorzeichen := positiv,
23// Ergebnis-Exponent := ceiling(e/2),
24// Ergebnis-Mantisse:
25//   Bilde aus [1,m22,...,m0,(26 Nullbits)] bei geradem e,
26//         aus [0,1,m22,...,m0,(25 Nullbits)] bei ungeradem e
27//   die Ganzzahl-Wurzel, eine 25-Bit-Zahl mit einer führenden 1.
28//   Runde das letzte Bit weg:
29//     Bit 0 = 0 -> abrunden,
30//     Bit 0 = 1 und Wurzel exakt -> round-to-even,
31//     Bit 0 = 1 und Rest >0 -> aufrunden.
32//   Dabei um ein Bit nach rechts schieben.
33//   Bei Aufrundung auf 2^24 (rounding overflow) Mantisse um 1 Bit nach rechts
34//     schieben und Exponent incrementieren.
35      // x entpacken:
36      var sintL exp;
37      var uint32 mant;
38      FF_decode(x, { return x; }, ,exp=,mant=);
39      // Um die 64-Bit-Ganzzahl-Wurzel ausnutzen zu können, fügen wir beim
40      // Radikanden 39 bzw. 40 statt 25 bzw. 26 Nullbits an.
41      if (exp & bit(0))
42        // e ungerade
43        { mant = mant << (31-(FF_mant_len+1)); exp = exp+1; }
44        else
45        // e gerade
46        { mant = mant << (32-(FF_mant_len+1)); }
47      exp = exp >> 1; // exp := exp/2
48      var bool exactp;
49      isqrt_64_32(mant,0, mant=,exactp=); // mant := isqrt(mant*2^32), eine 32-Bit-Zahl
50      // Die hinteren 31-FF_mant_len Bits wegrunden:
51      if ( ((mant & bit(30-FF_mant_len)) ==0) // Bit 7 =0 -> abrunden
52           || ( ((mant & (bit(30-FF_mant_len)-1)) ==0) // Bit 7 =1 und Bits 6..0 >0 -> aufrunden
53                && exactp                   // Bit 7 =1 und Bits 6..0 =0, aber Rest -> aufrunden
54                // round-to-even, je nach Bit 8 :
55                && ((mant & bit(31-FF_mant_len)) ==0)
56         )    )
57        // abrunden
58        { mant = mant >> (31-FF_mant_len); }
59        else
60        // aufrunden
61        { mant = mant >> (31-FF_mant_len);
62          mant += 1;
63          if (mant >= bit(FF_mant_len+1)) // rounding overflow?
64            { mant = mant>>1; exp = exp+1; }
65        }
66      return encode_FF(0,exp,mant);
67}
68
69}  // namespace cln