7 #include "cln/lfloat.h"
13 #include "cl_LF_impl.h"
16 #include "cln/abort.h"
20 const cl_LF sqrt (const cl_LF& x)
23 // x = 0.0 -> Ergebnis 0.0
24 // Ergebnis-Vorzeichen := positiv,
25 // Ergebnis-Exponent := ceiling(e/2),
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
32 // Runde das letzte Digit weg:
33 // Bit 15 = 0 -> abrunden,
34 // Bit 15 = 1, Rest =0 und Wurzel exakt -> round-to-even,
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;
45 var uintC r_len = 2*len+2; // Länge des Radikanden
46 num_stack_alloc(r_len, r_MSDptr=,r_LSDptr=);
47 if ((uexp & bit(0)) == (LF_exp_mid & bit(0)))
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
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
61 // Compute ((uexp - LF_exp_mid + 1) >> 1) + LF_exp_mid without risking
63 uexp = ((uexp - ((LF_exp_mid - 1) & 1)) >> 1) - ((LF_exp_mid - 1) >> 1)
65 // Ergebnis allozieren:
66 var Lfloat y = allocate_lfloat(len,uexp,0);
67 var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
69 #ifndef CL_LF_PEDANTIC
70 if (len > 2900) // This is about 15% faster
71 { // Kehrwert der Wurzel errechnen:
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:
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)
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
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
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
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
119 && ((lspref(w.LSDptr,1) & bit(0)) ==0)
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
132 // Bit complexity (N := length(x)): O(M(N)).