// 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)
{
// 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);
// 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)
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
}
}
#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
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
} }
}
// Bit complexity (N := length(x)): O(M(N)).
+} // namespace cln