4 #include "base/cl_sysdep.h"
7 #include "base/digitseq/cl_2DS.h"
12 #include "base/digit/cl_2D.h"
13 #include "base/digitseq/cl_DS.h"
15 // Break-even point of the Newton iteration vs. standard div2adic.
17 const unsigned int recip2adic_threshold = 620;
19 // Use the old default values from CLN version <= 1.0.3 as a crude estimate.
20 const unsigned int recip2adic_threshold = 380;
25 void recip2adic (uintC len, const uintD* a_LSDptr, uintD* dest_LSDptr)
28 // If len < threshold, use regular 2-adic division.
29 // Else [Newton iteration] set n := ceiling(len/2),
30 // compute recursively b := recip2adic(a mod 2^(intDsize*n)),
31 // return 2*b-a*b^2 mod 2^(intDsize*2*n).
33 var uintL k = 0; // number of Newton steps
35 while (n >= recip2adic_threshold) {
40 var uintD* one_LSDptr;
41 num_stack_alloc(n,,one_LSDptr=);
42 lspref(one_LSDptr,0) = 1;
43 clear_loop_lsp(one_LSDptr lspop 1,n-1);
44 div2adic(n,one_LSDptr,a_LSDptr,dest_LSDptr);
48 var uintD* prod_LSDptr;
49 num_stack_alloc(len+1,,b2_LSDptr=);
50 num_stack_alloc(2*len,,prod_LSDptr=);
52 // n = ceiling(len/2^k)
53 // Compute n2 = ceiling(len/2^(k-1)),
54 // then n = ceiling(n2/2).
56 var uintC n2 = ((len-1)>>k)+1; // = 2*n or = 2*n-1
57 // Set b := 2*b-a*b^2 mod 2^(intDsize*n2)
58 cl_UDS_mul_square(dest_LSDptr,n,b2_LSDptr); // b^2
59 cl_UDS_mul(b2_LSDptr,n2,a_LSDptr,n2,prod_LSDptr); // a*b^2
60 clear_loop_lsp(dest_LSDptr lspop n,n2-n);
61 shift1left_loop_lsp(dest_LSDptr,n+1); // (n+1 instead of n2 is ok)
62 subfrom_loop_lsp(prod_LSDptr,dest_LSDptr,n2);
67 // Bit complexity (N := len): O(M(N)).