15 // We observe the following timings:
16 // Time for square root of a_len = 2*N by b_len = N digits,
17 // on a i486 33 MHz running Linux:
29 // -----> Newton faster for 1570 <= N <= 1790 and for N >= 2100.
30 // When in doubt, prefer to choose the standard algorithm.
31 static inline cl_boolean cl_recipsqrt_suitable (uintL n)
32 { return (cl_boolean)(n >= 2100); }
34 // Workaround gcc-2.7.0 bug on i386.
37 #if (__GNUC_MINOR__ == 7)
38 #define workaround_gcc_bug() *&b_stern = *&b_stern;
42 #ifndef workaround_gcc_bug
43 #define workaround_gcc_bug()
46 // Bildet zu einer Unsigned Digit sequence a die Wurzel
47 // (genauer: Gaußklammer aus Wurzel aus a).
48 // squarep = cl_UDS_sqrt(a_MSDptr,a_len,a_LSDptr, &b);
49 // > a_MSDptr/a_len/a_LSDptr: eine UDS
50 // < NUDS b: Gaußklammer der Wurzel aus a
51 // < squarep: cl_true falls a = b^2, cl_false falls b^2 < a < (b+1)^2.
53 // erst A normalisieren. A=0 --> B=0, fertig.
54 // Wähle n so, daß beta^(2n-2) <= A < beta^(2n).
55 // Wähle s (0<=s<16) so, daß beta^(2n)/4 <= A*2^(2s) < beta^(2n).
56 // Setze A:=A*2^(2s) und kopiere dabei A. Suche B=floor(sqrt(A)).
57 // Mache Platz für B=[0,b[n-1],...,b[0]], (mit einem Nulldigit Platz davor,
58 // da dort nicht B, sondern 2*B abgespeichert werden wird).
59 // Auf den Plätzen [a[2n-1],...,a[2n-2j]] wird die Differenz
60 // [a[2n-1],...,a[2n-2j]] - [b[n-1],...,b[n-j]] ^ 2 abgespeichert.
61 // Bestimme b[n-1] = floor(sqrt(a[2n-1]*beta+a[2n-2])) mit Heron/Newton:
62 // {x:=beta als vorheriger Anfangswert, dann:}
63 // x := floor((beta+a[2n-1])/2)
64 // wiederhole: d:=floor((a[2n-1]*beta+a[2n-2])/x).
65 // Falls d<beta (kein Überlauf) und d<x,
66 // setze x:=floor((x+d)/2), nochmals.
67 // b[n-1]:=x. In B um ein Bit nach links verschoben abspeichern.
68 // {Wegen a[2n-1]>=beta/4 ist b[n-1]>=beta/2.}
69 // Erniedrige [a[2n-1],a[2n-2]] um b[n-1]^2.
71 // {Hier [b[n-1],...,b[n-j]] = floor(sqrt(altes [a[2n-1],...,a[2n-2j]])),
72 // in [a[2n-1],...,a[2n-2j]] steht jetzt der Rest
73 // [a[2n-1],...,a[2n-2j]] - [b[n-1],...,b[n-j]]^2, er ist >=0 und
74 // und <= 2 * [b[n-1],...,b[n-j]], belegt daher höchstens j Digits und 1 Bit.
75 // Daher sind nur [a[2n-j],...,a[2n-2j]] von Belang.}
76 // Für j<n: Bestimme die nächste Ziffer:
77 // b* := min(beta-1,floor([a[2n-j],...,a[2n-2j-1]]/(2*[b[n-1],...,b[n-j]]))).
78 // und [a[2n-j],...,a[2n-2j-1]] :=
79 // [a[2n-j],...,a[2n-2j-1]] - b* * 2 * [b[n-1],...,b[n-j]] (>= 0).
81 // b* := min(beta-1,floor([a[2n-j],a[2n-j-1],a[2n-j-2]]/(2*b[n-1]))),
82 // [a[2n-j],...,a[2n-2j-1]] wie angegeben erniedigen.
83 // Solange die Differenz <0 ist, setze b* := b* - 1 und
84 // erhöhe [a[2n-j],...,a[2n-2j-1]] um 2 * [b[n-1],...,b[n-j]].
85 // Erniedrige [a[2n-j],...,a[2n-2j-2]] um b* ^ 2.
86 // Tritt dabei ein negativer Carry auf,
87 // so setze b* := b* - 1,
88 // setze b[n-j-1] := b* (im Speicher um 1 Bit nach links verschoben),
89 // erhöhe [a[2n-j],...,a[2n-2j-2]] um 2*[b[n-1],...,b[n-j-1]]+1.
90 // Sonst setze b[n-j-1] := b* (im Speicher um 1 Bit nach links verschoben).
93 // Falls [a[n],...,a[0]] = [0,...,0], ist die Wurzel exakt, sonst nicht.
94 // Ergebnis ist [b[n-1],...,b[0]] * 2^(-s), schiebe also im Speicher
95 // [b[n],...,b[0]] um s+1 Bits nach rechts.
96 // Das Ergebnis ist eine NUDS der Länge n.
97 cl_boolean cl_UDS_sqrt (const uintD* a_MSDptr, uintC a_len, const uintD* a_LSDptr, DS* b_)
100 while ((a_len>0) && (mspref(a_MSDptr,0)==0)) { msshrink(a_MSDptr); a_len--; }
101 if (a_len==0) // A=0 -> B := NUDS 0
102 { b_->LSDptr = b_->MSDptr; b_->len = 0; return cl_true; }
104 // n und s bestimmen:
105 var uintC n = ceiling(a_len,2); // a_len = 2n oder 2n-1, n>0.
107 { var uintD msd = mspref(a_MSDptr,0); // a[2n] bzw. a[2n-1]
110 while /* ((msd & (bit(intDsize-1)|bit(intDsize-2))) ==0) */
111 (((sintD)msd >= 0) && ((sintD)(msd<<1) >= 0))
112 { msd = msd<<2; s++; }
114 integerlengthD(msd, s = intDsize - ); s = s>>1;
117 // Noch ist s nur modulo intDsize/2 bestimmt.
118 // A um 2s Bits nach links verschoben kopieren:
119 var uintD* new_a_MSDptr;
120 { var uintD* new_a_LSDptr;
121 num_stack_alloc(2*(uintL)n,new_a_MSDptr=,new_a_LSDptr=); // 2n Digits Platz belegen
122 {var uintL shiftcount = 2*s;
123 if (!((a_len & bit(0)) ==0)) // a_len ungerade?
124 { s += intDsize/2; lsprefnext(new_a_LSDptr) = 0; } // ja -> ein Nulldigit einschieben
126 { copy_loop_lsp(a_LSDptr,new_a_LSDptr,a_len); }
128 { shiftleftcopy_loop_lsp(a_LSDptr,new_a_LSDptr,a_len,shiftcount); }
130 #define a_MSDptr new_a_MSDptr
131 // Nun ist A = a_MSDptr/2n/..
132 if (cl_recipsqrt_suitable(n))
133 { // C := 1/sqrt(A) und dann D := A*C näherungsweise errechnen.
134 // D evtl. korrigieren, liefert B.
139 var uintD* d2_MSDptr;
140 num_stack_alloc(n+2, c_MSDptr=,c_LSDptr=);
141 num_stack_alloc(2*n+3, d_MSDptr=,d_LSDptr=);
142 num_stack_alloc(2*n, d2_MSDptr=,);
144 cl_UDS_recipsqrt(a_MSDptr,2*n,c_MSDptr,n);
145 // 1 <= c <= 2, | 1/sqrt(a) - c | < 1/2*beta^-n.
146 cl_UDS_mul(a_MSDptr mspop (n+1),n+1,c_LSDptr,n+2,d_LSDptr);
147 // 1/4 <= d < 2, | sqrt(a) - d | < beta^-n.
148 if (mspref(d_MSDptr,0) > 0)
149 { dec_loop_lsp(d_MSDptr mspop (n+1),n+1);
150 if (mspref(d_MSDptr,0) > 0) cl_abort();
152 // D is our guess for B. Square to see how much we have to correct.
153 cl_UDS_mul_square(d_MSDptr mspop (1+n),n,d2_MSDptr mspop 2*n);
155 b_->LSDptr = copy_loop_msp(d_MSDptr mspop 1,b_->MSDptr,n);
157 // Store 2*D in place of D.
158 if (shift1left_loop_lsp(d_MSDptr mspop (1+n),n))
159 mspref(d_MSDptr,0) = 1;
160 // Compare D^2 against A.
161 if (subfrom_loop_lsp(d2_MSDptr mspop 2*n,a_MSDptr mspop 2*n,2*n))
162 // guessed too high, decrement D
163 { dec_loop_lsp(b_->LSDptr,n);
164 dec_loop_lsp(d_MSDptr mspop (1+n),1+n); // store 2*D+1
165 if (!addto_loop_lsp(d_MSDptr mspop (1+n),a_MSDptr mspop 2*n,1+n))
167 if (!inc_loop_lsp(a_MSDptr mspop (n-1),n-1))
170 else if (test_loop_msp(a_MSDptr,n-1))
171 // guessed way too low
173 else if (compare_loop_msp(a_MSDptr mspop (n-1),d_MSDptr,1+n) > 0)
174 // guessed too low, increment D
175 { inc_loop_lsp(b_->LSDptr,n);
176 mspref(d_MSDptr,n) |= bit(0); // store 2*D-1
177 subfrom_loop_lsp(d_MSDptr mspop (1+n),a_MSDptr mspop 2*n,1+n);
178 inc_loop_lsp(d_MSDptr mspop (1+n),1+n); // store 2*D
179 if (compare_loop_msp(a_MSDptr mspop (n-1),d_MSDptr,1+n) > 0)
185 // Schiebe b um s Bits nach rechts:
187 shiftright_loop_msp(b_->MSDptr,n,s);
188 // Teste, ob alle a[n],...,a[0]=0 sind:
189 if (test_loop_msp(a_MSDptr mspop (n-1),n+1))
192 return cl_true; // ja -> Wurzel exakt
194 // Platz für B belegen:
195 { var uintD* b_MSDptr = b_->MSDptr mspop -1; // ab hier n+1 Digits Platz
197 // B = [0,b[n-1],...,b[0]] = b_MSDptr/n+1/..
198 // Bestimmung von b[n-1]:
199 { var uintD a_msd = mspref(a_MSDptr,0); // a[2n-1]
200 var uintD a_2msd = mspref(a_MSDptr,1); // a[2n-2]
202 var uintDD a_msdd = highlowDD(a_msd,a_2msd); // a[2n-1]*beta+a[2n-2]
204 // Anfangswert: x := floor((beta + a[2n-1])/2)
205 var uintD x = floor(a_msd,2) | bit(intDsize-1);
206 loop // Heron-Iterationsschleife
208 // Dividiere d := floor((a[2n-1]*beta+a[2n-2])/x) :
209 if (a_msd>=x) break; // Überlauf -> d>=beta -> fertig
211 divuD(a_msdd,x, d=,);
213 divuD(a_msd,a_2msd,x, d=,);
215 if (d >= x) break; // d>=x -> fertig
216 // Nächste Iteration: x := floor((x+d)/2)
217 // (Da die Folge der x bekanntlich monoton fallend ist
218 // und bei b[n-1] >= beta/2 endet, muß x >= beta/2 werden,
221 x = (uintD)(floor((uintDD)x + (uintDD)d, 2));
223 x = floor((uintD)(x+d),2) | bit(intDsize-1);
226 // x = b[n-1] fertig berechnet.
228 // Quadrieren und von [a[2n-1],a[2n-2]] abziehen:
230 a_msdd -= muluD(x,x);
231 mspref(a_MSDptr,0) = highD(a_msdd); mspref(a_MSDptr,1) = lowD(a_msdd);
235 muluD(x,x, x2hi=,x2lo=);
236 mspref(a_MSDptr,0) = a_msd - x2hi;
237 if (a_2msd < x2lo) { mspref(a_MSDptr,0) -= 1; }
238 mspref(a_MSDptr,1) = a_2msd - x2lo;
241 mspref(b_MSDptr,0) = 1; mspref(b_MSDptr,1) = x<<1; // b[n-1] ablegen
244 var uintD* a_mptr = a_MSDptr mspop 0;
245 var uintD* a_lptr = a_MSDptr mspop 2;
246 var uintD* b_ptr = b_MSDptr mspop 2;
247 // Wurzel-Hauptschleife
248 until (++j == n) // j=1,...,n
249 { // b_MSDptr = Pointer auf b[n], b_ptr = Pointer hinter b[n-j].
250 // a_mptr = Pointer auf a[2n-j], a_lptr = Pointer hinter a[2n-2j].
253 { var uintD a_1d = mspref(a_mptr,0); // a[2n-j], =0 oder =1
254 var uintD a_2d = mspref(a_mptr,1); // a[2n-j-1]
255 var uintD a_3d = mspref(a_mptr,2); // a[2n-j-2]
256 // a[2n-j]*beta^2+a[2n-j-1]*beta+a[2n-j-2] durch 2 dividieren,
257 // dann durch b_msd = b[n-1] dividieren:
259 var uintDD a_123dd = highlowDD(a_2d,a_3d);
260 a_123dd = a_123dd>>1; if (!(a_1d==0)) { a_123dd |= bit(2*intDsize-1); }
261 if (highD(a_123dd) >= b_msd)
262 { b_stern = bitm(intDsize)-1; } // bei Überlauf: beta-1
264 { divuD(a_123dd,b_msd, b_stern=,); }
266 a_3d = a_3d>>1; if (!((a_2d & bit(0)) ==0)) { a_3d |= bit(intDsize-1); }
267 a_2d = a_2d>>1; if (!(a_1d==0)) { a_2d |= bit(intDsize-1); }
269 { b_stern = bitm(intDsize)-1; } // bei Überlauf: beta-1
271 { divuD(a_2d,a_3d,b_msd, b_stern=,); }
274 // b_stern = b* in der ersten Schätzung.
275 a_lptr = a_lptr mspop 1; // Pointer hinter a[2n-2j-1]
276 // Subtraktion [a[2n-j],...,a[2n-2j-1]] -= b* * [b[n],b[n-1],...,b[n-j]] :
277 { var uintD carry = mulusub_loop_lsp(b_stern,b_ptr,a_lptr,j+1);
278 if (mspref(a_mptr,0) >= carry)
279 { mspref(a_mptr,0) -= carry; }
281 { mspref(a_mptr,0) -= carry; // a[2n-j] wird <0
282 // negativer Übertrag -> b* nach unten korrigieren:
284 { b_stern = b_stern-1; // b* := b* - 1
285 // erhöhe [a[2n-j],...,a[2n-2j-1]] um [b[n],...,b[n-j]]:
286 if (!(( addto_loop_lsp(b_ptr,a_lptr,j+1) ==0)))
287 if ((mspref(a_mptr,0) += 1) ==0) // Übertrag zu a[2n-j]
288 break; // macht a[2n-j] wieder >=0 -> Subtraktionsergebnis >=0
290 // b_stern = b* in der zweiten Schätzung.
291 a_mptr = a_mptr mspop 1; // Pointer auf a[2n-j-1]
292 a_lptr = a_lptr mspop 1; // Pointer hinter a[2n-2j-2]
293 // Ziehe b* ^ 2 von [a[2n-j],...,a[2n-2j-2]] ab:
295 { var uintDD b_stern_2 = muluD(b_stern,b_stern);
296 var uintDD a_12dd = highlowDD(lspref(a_lptr,1),lspref(a_lptr,0)); // a[2n-2j-1]*beta+a[2n-2j-2]
297 var uintDD a_12dd_new = a_12dd - b_stern_2;
298 lspref(a_lptr,1) = highD(a_12dd_new); lspref(a_lptr,0) = lowD(a_12dd_new);
299 if (a_12dd >= b_stern_2) goto b_stern_ok;
302 { var uintD b_stern_2_hi;
303 var uintD b_stern_2_lo;
304 muluD(b_stern,b_stern, b_stern_2_hi=,b_stern_2_lo=);
305 {var uintD a_1d = lspref(a_lptr,1); // a[2n-2j-1]
306 var uintD a_2d = lspref(a_lptr,0); // a[2n-2j-2]
307 var uintD a_1d_new = a_1d - b_stern_2_hi;
308 var uintD a_2d_new = a_2d - b_stern_2_lo;
309 if (a_2d < b_stern_2_lo) { a_1d_new -= 1; }
310 lspref(a_lptr,1) = a_1d_new; lspref(a_lptr,0) = a_2d_new;
311 if ((a_1d > b_stern_2_hi)
312 || ((a_1d == b_stern_2_hi) && (a_2d >= b_stern_2_lo))
318 { // muß noch [a[2n-j],...,a[2n-2j]] um 1 erniedrigen:
319 if ( dec_loop_lsp(a_lptr lspop 2,j+1) ==0) goto b_stern_ok;
320 // Subtraktion von b*^2 lieferte negativen Carry
321 b_stern = b_stern-1; // b* := b* - 1
322 workaround_gcc_bug();
323 // erhöhe [a[2n-j-1],...,a[2n-2j-2]] um [b[n],...,b[n-j],0] + 2 * b* + 1
324 if ((sintD)b_stern < 0) { mspref(b_ptr,-1) |= bit(0); } // höchstes Bit von b* in b[n-j] ablegen
325 mspref(b_ptr,0) = (uintD)(b_stern<<1)+1; // niedrige Bits von b* und eine 1 als b[n-j-1] ablegen
326 addto_loop_lsp(b_ptr mspop 1,a_lptr,j+2);
327 // (a[2n-j] wird nicht mehr gebraucht.)
328 mspref(b_ptr,0) -= 1; // niedrige Bits von b* in b[n-j-1] ablegen
329 b_ptr = b_ptr mspop 1;
333 { // b* als b[n-j-1] ablegen:
334 if ((sintD)b_stern < 0) { mspref(b_ptr,-1) |= bit(0); } // höchstes Bit von b* in b[n-j] ablegen
335 mspref(b_ptr,0) = (uintD)(b_stern<<1); // niedrige Bits von b* als b[n-j-1] ablegen
336 b_ptr = b_ptr mspop 1;
339 // b_MSDptr = Pointer auf b[n], b_ptr = Pointer hinter b[0].
340 // a_mptr = Pointer auf a[n].
341 // Schiebe [b[n],...,b[0]] um s+1 Bits nach rechts:
345 { shiftright_loop_msp(b_MSDptr,n+1,s+1); msshrink(b_MSDptr); }
346 // b = b_MSDptr/n/b_ptr ist fertig, eine NUDS.
347 b_->MSDptr = b_MSDptr; b_->len = n; b_->LSDptr = b_ptr;
348 // Teste, ob alle a[n],...,a[0]=0 sind:
349 if (test_loop_msp(a_mptr,n+1))
352 { return cl_true; } // ja -> Wurzel exakt
355 // Bit complexity (N := a_len): O(M(N)).