PageRenderTime 63ms CodeModel.GetById 21ms RepoModel.GetById 1ms app.codeStats 0ms

/cln-1.3.2/src/base/digitseq/cl_DS_sqrt.cc

#
C++ | 364 lines | 221 code | 9 blank | 134 comment | 46 complexity | d193bd2d3f61caa44d4e6a9e567e18c9 MD5 | raw file
Possible License(s): GPL-2.0
  1. // cl_UDS_sqrt().
  2. // General includes.
  3. #include "base/cl_sysdep.h"
  4. // Specification.
  5. #include "base/digitseq/cl_DS.h"
  6. // Implementation.
  7. #include "base/cl_low.h"
  8. #include "cln/exception.h"
  9. namespace cln {
  10. // We observe the following timings:
  11. // Time for square root of a_len = 2*N by b_len = N digits,
  12. // OS: Linux 2.2, intDsize==32, OS: TRU64/4.0, intDsize==64,
  13. // Machine: P-III/450MHz Machine: EV5/300MHz:
  14. // N standard Newton standard Newton
  15. // 30 0.00002 0.00009 0.00011 0.00027
  16. // 100 0.00012 0.00052 0.00057 0.0017
  17. // 300 0.00087 0.0031 0.0037 0.0091
  18. // 1000 0.0089 0.020 0.037 0.069
  19. // 3000 0.087 0.11 <-(~3200) 0.30 0.28 <- (~2750)
  20. // 10000 1.27 0.55 3.5 1.3
  21. // 30000 12.7 1.35 31.1 3.4
  22. // Newton faster for 3200<N Newton faster for 2750<N
  23. // When in doubt, prefer to choose the standard algorithm.
  24. #if CL_USE_GMP
  25. static inline bool cl_recipsqrt_suitable (uintC n)
  26. { return n >= 3200; }
  27. #else
  28. // Use the old default values from CLN version <= 1.0.3 as a crude estimate.
  29. // Time for square root of a_len = 2*N by b_len = N digits,
  30. // on a i486 33 MHz running Linux:
  31. // N standard Newton
  32. // 10 0.00022 0.00132
  33. // 25 0.00082 0.0047
  34. // 50 0.0026 0.0130
  35. // 100 0.0095 0.038
  36. // 250 0.057 0.154
  37. // 500 0.22 0.46
  38. // 1000 0.90 1.39
  39. // 2500 6.0 4.6
  40. // 5000 24.1 10.7
  41. // 10000 98 23.2
  42. // -----> Newton faster for 1570 <= N <= 1790 and for N >= 2100.
  43. static inline bool cl_recipsqrt_suitable (uintC n)
  44. { return n >= 2100; }
  45. #endif
  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: true falls a = b^2, false falls b^2 < a < (b+1)^2.
  52. // Methode:
  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.
  70. // Für j=1,...,n:
  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).
  80. // Im einzelnen:
  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).
  91. // Nächstes j.
  92. // Für j=n:
  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. bool cl_UDS_sqrt (const uintD* a_MSDptr, uintC a_len, const uintD* a_LSDptr, DS* b_)
  98. {
  99. // A normalisieren:
  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 true; }
  103. CL_ALLOCA_STACK;
  104. // n und s bestimmen:
  105. var uintC n = ceiling(a_len,2); // a_len = 2n oder 2n-1, n>0.
  106. var uintL s;
  107. { var uintD msd = mspref(a_MSDptr,0); // a[2n] bzw. a[2n-1]
  108. #if 0
  109. s = 0;
  110. while /* ((msd & (bit(intDsize-1)|bit(intDsize-2))) ==0) */
  111. (((sintD)msd >= 0) && ((sintD)(msd<<1) >= 0))
  112. { msd = msd<<2; s++; }
  113. #else
  114. integerlengthD(msd, s = intDsize - ); s = s>>1;
  115. #endif
  116. }
  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*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
  125. if (shiftcount==0)
  126. { copy_loop_lsp(a_LSDptr,new_a_LSDptr,a_len); }
  127. else
  128. { shiftleftcopy_loop_lsp(a_LSDptr,new_a_LSDptr,a_len,shiftcount); }
  129. }}
  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.
  135. var uintD* c_MSDptr;
  136. var uintD* c_LSDptr;
  137. var uintD* d_MSDptr;
  138. var uintD* d_LSDptr;
  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=,);
  143. // 1/4 <= a < 1.
  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) throw runtime_exception();
  151. }
  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);
  154. // Store D.
  155. b_->LSDptr = copy_loop_msp(d_MSDptr mspop 1,b_->MSDptr,n);
  156. b_->len = 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))
  166. throw runtime_exception();
  167. if (!inc_loop_lsp(a_MSDptr mspop (n-1),n-1))
  168. throw runtime_exception();
  169. }
  170. else if (test_loop_msp(a_MSDptr,n-1))
  171. // guessed way too low
  172. throw runtime_exception();
  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)
  180. throw runtime_exception();
  181. }
  182. else
  183. // guessed ok
  184. {}
  185. // Schiebe b um s Bits nach rechts:
  186. if (s > 0)
  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))
  190. return false;
  191. else
  192. return true; // ja -> Wurzel exakt
  193. }
  194. // Platz für B belegen:
  195. { var uintD* b_MSDptr = b_->MSDptr mspop -1; // ab hier n+1 Digits Platz
  196. var uintD b_msd;
  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]
  201. #if HAVE_DD
  202. var uintDD a_msdd = highlowDD(a_msd,a_2msd); // a[2n-1]*beta+a[2n-2]
  203. #endif
  204. // Anfangswert: x := floor((beta + a[2n-1])/2)
  205. var uintD x = floor(a_msd,2) | bit(intDsize-1);
  206. loop // Heron-Iterationsschleife
  207. { var uintD d;
  208. // Dividiere d := floor((a[2n-1]*beta+a[2n-2])/x) :
  209. if (a_msd>=x) break; // Überlauf -> d>=beta -> fertig
  210. #if HAVE_DD
  211. divuD(a_msdd,x, d=,);
  212. #else
  213. divuD(a_msd,a_2msd,x, d=,);
  214. #endif
  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,
  219. // d.h. x+d>=beta.)
  220. #if HAVE_DD
  221. x = (uintD)(floor((uintDD)x + (uintDD)d, 2));
  222. #else
  223. x = floor((uintD)(x+d),2) | bit(intDsize-1);
  224. #endif
  225. }
  226. // x = b[n-1] fertig berechnet.
  227. b_msd = x;
  228. // Quadrieren und von [a[2n-1],a[2n-2]] abziehen:
  229. #if HAVE_DD
  230. a_msdd -= muluD(x,x);
  231. mspref(a_MSDptr,0) = highD(a_msdd); mspref(a_MSDptr,1) = lowD(a_msdd);
  232. #else
  233. {var uintD x2hi;
  234. var uintD x2lo;
  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;
  239. }
  240. #endif
  241. mspref(b_MSDptr,0) = 1; mspref(b_MSDptr,1) = x<<1; // b[n-1] ablegen
  242. }
  243. {var uintC j = 0;
  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].
  251. // Bestimme b* :
  252. var uintD b_stern;
  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:
  258. #if HAVE_DD
  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
  263. else
  264. { divuD(a_123dd,b_msd, b_stern=,); }
  265. #else
  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); }
  268. if (a_2d >= b_msd)
  269. { b_stern = bitm(intDsize)-1; } // bei Überlauf: beta-1
  270. else
  271. { divuD(a_2d,a_3d,b_msd, b_stern=,); }
  272. #endif
  273. }
  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; }
  280. else
  281. { mspref(a_mptr,0) -= carry; // a[2n-j] wird <0
  282. // negativer Übertrag -> b* nach unten korrigieren:
  283. loop
  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
  289. } } }
  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:
  294. #if HAVE_DD
  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;
  300. }
  301. #else
  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))
  313. )
  314. goto b_stern_ok;
  315. }}
  316. #endif
  317. if (TRUE)
  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. // erhöhe [a[2n-j-1],...,a[2n-2j-2]] um [b[n],...,b[n-j],0] + 2 * b* + 1
  323. if ((sintD)b_stern < 0) { mspref(b_ptr,-1) |= bit(0); } // höchstes Bit von b* in b[n-j] ablegen
  324. mspref(b_ptr,0) = (uintD)(b_stern<<1)+1; // niedrige Bits von b* und eine 1 als b[n-j-1] ablegen
  325. addto_loop_lsp(b_ptr mspop 1,a_lptr,j+2);
  326. // (a[2n-j] wird nicht mehr gebraucht.)
  327. mspref(b_ptr,0) -= 1; // niedrige Bits von b* in b[n-j-1] ablegen
  328. b_ptr = b_ptr mspop 1;
  329. }
  330. else
  331. b_stern_ok:
  332. { // b* als b[n-j-1] ablegen:
  333. if ((sintD)b_stern < 0) { mspref(b_ptr,-1) |= bit(0); } // höchstes Bit von b* in b[n-j] ablegen
  334. mspref(b_ptr,0) = (uintD)(b_stern<<1); // niedrige Bits von b* als b[n-j-1] ablegen
  335. b_ptr = b_ptr mspop 1;
  336. }
  337. }
  338. // b_MSDptr = Pointer auf b[n], b_ptr = Pointer hinter b[0].
  339. // a_mptr = Pointer auf a[n].
  340. // Schiebe [b[n],...,b[0]] um s+1 Bits nach rechts:
  341. if (s == intDsize-1)
  342. { lsshrink(b_ptr); }
  343. else
  344. { shiftright_loop_msp(b_MSDptr,n+1,s+1); msshrink(b_MSDptr); }
  345. // b = b_MSDptr/n/b_ptr ist fertig, eine NUDS.
  346. b_->MSDptr = b_MSDptr; b_->len = n; b_->LSDptr = b_ptr;
  347. // Teste, ob alle a[n],...,a[0]=0 sind:
  348. if (test_loop_msp(a_mptr,n+1))
  349. { return false; }
  350. else
  351. { return true; } // ja -> Wurzel exakt
  352. }}
  353. }
  354. // Bit complexity (N := a_len): O(M(N)).
  355. } // namespace cln