src/float/lfloat/algebraic/cl_LF_sqrt.cc

   1 // sqrt().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/lfloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_LF.h"
  13 #include "cl_LF_impl.h"
  14 #include "cl_F.h"
  15 #include "cl_DS.h"
  16 #include "cln/abort.h"
  17
  18 namespace cln {
  19
  20 const cl_LF sqrt (const cl_LF& x)
  21 {
  22 // Methode:
  23 // x = 0.0 -> Ergebnis 0.0
  24 // Ergebnis-Vorzeichen := positiv,
  25 // Ergebnis-Exponent := ceiling(e/2),
  26 // Ergebnis-Mantisse:
  27 //   Erweitere die Mantisse (n Digits) um n+2 Nulldigits nach hinten.
  28 //   Bei ungeradem e schiebe dies (oder nur die ersten n+1 Digits davon)
  29 //     um 1 Bit nach rechts.
  30 //   Bilde daraus die Ganzzahl-Wurzel, eine n+1-Digit-Zahl mit einer
  31 //     fhrenden 1.
  32 //   Runde das letzte Digit weg:
  33 //     Bit 15 = 0 -> abrunden,
  34 //     Bit 15 = 1, Rest =0 und Wurzel exakt -> round-to-even,
  35 //     sonst aufrunden.
  36 //   Bei rounding overflow Mantisse um 1 Bit nach rechts schieben
  37 //     und Exponent incrementieren.
  38       var uintL uexp = TheLfloat(x)->expo;
  39       if (uexp==0) { return x; } // x=0.0 -> 0.0 als Ergebnis
  40       var uintC len = TheLfloat(x)->len;
  41       // Radikanden bilden:
  42       CL_ALLOCA_STACK;
  43       var uintD* r_MSDptr;
  44       var uintD* r_LSDptr;
  45       var uintL r_len = 2*(uintL)len+2; // L�ge des Radikanden
  46       num_stack_alloc(r_len, r_MSDptr=,r_LSDptr=);
  47       if ((uexp & bit(0)) == (LF_exp_mid & bit(0)))
  48         // Exponent gerade
  49         {var uintD* ptr =
  50            copy_loop_msp(arrayMSDptr(TheLfloat(x)->data,len),r_MSDptr,len); // n Digits kopieren
  51          clear_loop_msp(ptr,len+2); // n+2 Nulldigits anh�gen
  52         }
  53         else
  54         // Exponent ungerade
  55         {var uintD carry_rechts = // n Digits kopieren und um 1 Bit rechts shiften
  56            shiftrightcopy_loop_msp(arrayMSDptr(TheLfloat(x)->data,len),r_MSDptr,len,1,0);
  57          var uintD* ptr = r_MSDptr mspop len;
  58          msprefnext(ptr) = carry_rechts; // �ertrag und
  59          clear_loop_msp(ptr,len+1); // n+1 Nulldigits anh�gen
  60         }
  61       // Compute ((uexp - LF_exp_mid + 1) >> 1) + LF_exp_mid without risking
  62       // uintL overflow.
  63       uexp = ((uexp - ((LF_exp_mid - 1) & 1)) >> 1) - ((LF_exp_mid - 1) >> 1)
  64              + LF_exp_mid;
  65       // Ergebnis allozieren:
  66       var Lfloat y = allocate_lfloat(len,uexp,0);
  67       var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
  68       // Wurzel ziehen:
  69 #ifndef CL_LF_PEDANTIC
  70       if (len > 2900) // This is about 15% faster
  71         { // Kehrwert der Wurzel errechnen:
  72           var uintD* s_MSDptr;
  73           var uintD* s_LSDptr;
  74           num_stack_alloc(len+2, s_MSDptr=,s_LSDptr=);
  75           cl_UDS_recipsqrt(r_MSDptr,r_len, s_MSDptr,len);
  76           // Mit dem Radikanden multiplizieren:
  77           var uintD* p_MSDptr;
  78           var uintD* p_LSDptr;
  79           num_stack_alloc(2*len+3, p_MSDptr=,p_LSDptr=);
  80           cl_UDS_mul(r_MSDptr mspop (len+1),len+1,s_LSDptr,len+2,p_LSDptr);
  81           // Ablegen und runden:
  82           copy_loop_msp(p_MSDptr mspop 1,y_mantMSDptr,len); // NUDS nach y kopieren
  83           if (mspref(p_MSDptr,0) == 0)
  84             { if ( ((sintD)mspref(p_MSDptr,len+1) >= 0) // n�hstes Bit =0 -> abrunden
  85                    || ( ((mspref(p_MSDptr,len+1) & ((uintD)bit(intDsize-1)-1)) ==0) // =1 und weitere Bits >0 -> aufrunden
  86                         && !test_loop_msp(p_MSDptr mspop (len+2),len+1)
  87                         // round-to-even (etwas witzlos, da eh alles ungenau ist)
  88                         && ((mspref(p_MSDptr,len) & bit(0)) ==0)
  89                  )    )
  90                 // abrunden
  91                 {}
  92                 else
  93                 // aufrunden
  94                 { if ( inc_loop_lsp(y_mantMSDptr mspop len,len) )
  95                     // �ertrag durchs Aufrunden
  96                     { mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
  97                       (TheLfloat(y)->expo)++; // Exponenten incrementieren
  98                 }   }
  99             }
 100             else
 101             // �ertrag durch Rundungsfehler
 102             { if (test_loop_msp(y_mantMSDptr,len)) cl_abort();
 103               mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
 104               (TheLfloat(y)->expo)++; // Exponenten incrementieren
 105             }
 106           return y;
 107         }
 108 #endif
 109       var DS w;
 110       var cl_boolean exactp;
 111       UDS_sqrt(r_MSDptr,r_len,r_LSDptr, &w, exactp=);
 112       // w ist die Ganzzahl-Wurzel, eine n+1-Digit-Zahl.
 113       copy_loop_msp(w.MSDptr,y_mantMSDptr,len); // NUDS nach y kopieren
 114       // Runden:
 115       if ( ((sintD)lspref(w.LSDptr,0) >= 0) // n�hstes Bit =0 -> abrunden
 116            || ( ((lspref(w.LSDptr,0) & ((uintD)bit(intDsize-1)-1)) ==0) // =1 und weitere Bits >0 oder Rest >0 -> aufrunden
 117                 && exactp
 118                 // round-to-even
 119                 && ((lspref(w.LSDptr,1) & bit(0)) ==0)
 120          )    )
 121         // abrunden
 122         {}
 123         else
 124         // aufrunden
 125         { if ( inc_loop_lsp(y_mantMSDptr mspop len,len) )
 126             // �ertrag durchs Aufrunden
 127             { mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
 128               (TheLfloat(y)->expo)++; // Exponenten incrementieren
 129         }   }
 130       return y;
 131 }
 132 // Bit complexity (N := length(x)): O(M(N)).
 133
 134 }  // namespace cln