src/float/sfloat/algebraic/cl_SF_sqrt.cc

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