PageRenderTime 45ms CodeModel.GetById 18ms RepoModel.GetById 1ms app.codeStats 0ms

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

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