Use paths relative the `src' directory in the #include directives.

[cln.git] / src / float / lfloat / algebraic / cl_LF_sqrt.cc

diff --git a/src/float/lfloat/algebraic/cl_LF_sqrt.cc b/src/float/lfloat/algebraic/cl_LF_sqrt.cc
index 9679a4076cd2deed7358663186594effa4f24646..b58f598c94c53b81c3e7487d1dc439c17c62abf5 100644 (file)
--- a/src/float/lfloat/algebraic/cl_LF_sqrt.cc
+++ b/src/float/lfloat/algebraic/cl_LF_sqrt.cc
@@ -1,19 +1,21 @@
 // sqrt().
 
 // General includes.
-#include "cl_sysdep.h"
+#include "base/cl_sysdep.h"
 
 // Specification.
-#include "cl_lfloat.h"
+#include "cln/lfloat.h"
 
 
 // Implementation.
 
-#include "cl_LF.h"
-#include "cl_LF_impl.h"
-#include "cl_F.h"
-#include "cl_DS.h"
-#include "cl_abort.h"
+#include "float/lfloat/cl_LF.h"
+#include "float/lfloat/cl_LF_impl.h"
+#include "float/cl_F.h"
+#include "base/digitseq/cl_DS.h"
+#include "cln/exception.h"
+
+namespace cln {
 
 const cl_LF sqrt (const cl_LF& x)
 {
@@ -26,45 +28,46 @@ const cl_LF sqrt (const cl_LF& x)
 //   Bei ungeradem e schiebe dies (oder nur die ersten n+1 Digits davon)
 //     um 1 Bit nach rechts.
 //   Bilde daraus die Ganzzahl-Wurzel, eine n+1-Digit-Zahl mit einer
-//     führenden 1.
+//     fhrenden 1.
 //   Runde das letzte Digit weg:
 //     Bit 15 = 0 -> abrunden,
 //     Bit 15 = 1, Rest =0 und Wurzel exakt -> round-to-even,
 //     sonst aufrunden.
 //   Bei rounding overflow Mantisse um 1 Bit nach rechts schieben
 //     und Exponent incrementieren.
-      var uintL uexp = TheLfloat(x)->expo;
+      var uintE uexp = TheLfloat(x)->expo;
       if (uexp==0) { return x; } // x=0.0 -> 0.0 als Ergebnis
       var uintC len = TheLfloat(x)->len;
       // Radikanden bilden:
       CL_ALLOCA_STACK;
       var uintD* r_MSDptr;
       var uintD* r_LSDptr;
-      var uintL r_len = 2*(uintL)len+2; // Länge des Radikanden
+      var uintC r_len = 2*len+2; // Länge des Radikanden
       num_stack_alloc(r_len, r_MSDptr=,r_LSDptr=);
-      uexp = uexp - LF_exp_mid + 1;
-      if (uexp & bit(0))
+      if ((uexp & bit(0)) == (LF_exp_mid & bit(0)))
         // Exponent gerade
         {var uintD* ptr =
            copy_loop_msp(arrayMSDptr(TheLfloat(x)->data,len),r_MSDptr,len); // n Digits kopieren
-         clear_loop_msp(ptr,len+2); // n+2 Nulldigits anhängen
+         clear_loop_msp(ptr,len+2); // n+2 Nulldigits anh�gen
         }
         else
         // Exponent ungerade
         {var uintD carry_rechts = // n Digits kopieren und um 1 Bit rechts shiften
            shiftrightcopy_loop_msp(arrayMSDptr(TheLfloat(x)->data,len),r_MSDptr,len,1,0);
          var uintD* ptr = r_MSDptr mspop len;
-         msprefnext(ptr) = carry_rechts; // Übertrag und
-         clear_loop_msp(ptr,len+1); // n+1 Nulldigits anhängen
+         msprefnext(ptr) = carry_rechts; // �ertrag und
+         clear_loop_msp(ptr,len+1); // n+1 Nulldigits anh�gen
         }
-      uexp = (sintL)((sintL)uexp >> 1); // Exponent halbieren
-      uexp = uexp + LF_exp_mid;
+      // Compute ((uexp - LF_exp_mid + 1) >> 1) + LF_exp_mid without risking
+      // uintE overflow.
+      uexp = ((uexp - ((LF_exp_mid - 1) & 1)) >> 1) - ((LF_exp_mid - 1) >> 1)
+             + LF_exp_mid;
       // Ergebnis allozieren:
       var Lfloat y = allocate_lfloat(len,uexp,0);
       var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
       // Wurzel ziehen:
 #ifndef CL_LF_PEDANTIC
-      if (len > 1900) // Das ist etwa 10% bis 20% schneller (im Mittel 15%).
+      if (len > 2900) // This is about 15% faster
         { // Kehrwert der Wurzel errechnen:
           var uintD* s_MSDptr;
           var uintD* s_LSDptr;
@@ -78,7 +81,7 @@ const cl_LF sqrt (const cl_LF& x)
           // Ablegen und runden:
           copy_loop_msp(p_MSDptr mspop 1,y_mantMSDptr,len); // NUDS nach y kopieren
           if (mspref(p_MSDptr,0) == 0)
-            { if ( ((sintD)mspref(p_MSDptr,len+1) >= 0) // nächstes Bit =0 -> abrunden
+            { if ( ((sintD)mspref(p_MSDptr,len+1) >= 0) // n�hstes Bit =0 -> abrunden
                    || ( ((mspref(p_MSDptr,len+1) & ((uintD)bit(intDsize-1)-1)) ==0) // =1 und weitere Bits >0 -> aufrunden
                         && !test_loop_msp(p_MSDptr mspop (len+2),len+1)
                         // round-to-even (etwas witzlos, da eh alles ungenau ist)
@@ -89,14 +92,14 @@ const cl_LF sqrt (const cl_LF& x)
                 else
                 // aufrunden
                 { if ( inc_loop_lsp(y_mantMSDptr mspop len,len) )
-                    // Übertrag durchs Aufrunden
+                    // �ertrag durchs Aufrunden
                     { mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
                       (TheLfloat(y)->expo)++; // Exponenten incrementieren
                 }   }
             }
             else
-            // Übertrag durch Rundungsfehler
-            { if (test_loop_msp(y_mantMSDptr,len)) cl_abort();
+            // �ertrag durch Rundungsfehler
+            { if (test_loop_msp(y_mantMSDptr,len)) throw runtime_exception();
               mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
               (TheLfloat(y)->expo)++; // Exponenten incrementieren
             }
@@ -104,12 +107,12 @@ const cl_LF sqrt (const cl_LF& x)
         }
 #endif
       var DS w;
-      var cl_boolean exactp;
+      var bool exactp;
       UDS_sqrt(r_MSDptr,r_len,r_LSDptr, &w, exactp=);
       // w ist die Ganzzahl-Wurzel, eine n+1-Digit-Zahl.
       copy_loop_msp(w.MSDptr,y_mantMSDptr,len); // NUDS nach y kopieren
       // Runden:
-      if ( ((sintD)lspref(w.LSDptr,0) >= 0) // nächstes Bit =0 -> abrunden
+      if ( ((sintD)lspref(w.LSDptr,0) >= 0) // n�hstes Bit =0 -> abrunden
            || ( ((lspref(w.LSDptr,0) & ((uintD)bit(intDsize-1)-1)) ==0) // =1 und weitere Bits >0 oder Rest >0 -> aufrunden
                 && exactp
                 // round-to-even
@@ -120,7 +123,7 @@ const cl_LF sqrt (const cl_LF& x)
         else
         // aufrunden
         { if ( inc_loop_lsp(y_mantMSDptr mspop len,len) )
-            // Übertrag durchs Aufrunden
+            // �ertrag durchs Aufrunden
             { mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
               (TheLfloat(y)->expo)++; // Exponenten incrementieren
         }   }
@@ -128,3 +131,4 @@ const cl_LF sqrt (const cl_LF& x)
 }
 // Bit complexity (N := length(x)): O(M(N)).
 
+}  // namespace cln