15 // Break-even point of the Newton iteration vs. standard div2adic.
16 const unsigned int recip2adic_threshold = 380;
18 void recip2adic (uintC len, const uintD* a_LSDptr, uintD* dest_LSDptr)
21 // If len < threshold, use regular 2-adic division.
22 // Else [Newton iteration] set n := ceiling(len/2),
23 // compute recursively b := recip2adic(a mod 2^(intDsize*n)),
24 // return 2*b-a*b^2 mod 2^(intDsize*2*n).
26 var uintL k = 0; // number of Newton steps
28 while (n >= recip2adic_threshold) {
33 var uintD* one_LSDptr;
34 num_stack_alloc(n,,one_LSDptr=);
35 lspref(one_LSDptr,0) = 1;
36 clear_loop_lsp(one_LSDptr lspop 1,n-1);
37 div2adic(n,one_LSDptr,a_LSDptr,dest_LSDptr);
41 var uintD* prod_LSDptr;
42 num_stack_alloc(len+1,,b2_LSDptr=);
43 num_stack_alloc(2*(uintL)len,,prod_LSDptr=);
45 // n = ceiling(len/2^k)
46 // Compute n2 = ceiling(len/2^(k-1)),
47 // then n = ceiling(n2/2).
49 var uintL n2 = ((len-1)>>k)+1; // = 2*n or = 2*n-1
50 // Set b := 2*b-a*b^2 mod 2^(intDsize*n2)
51 cl_UDS_mul_square(dest_LSDptr,n,b2_LSDptr); // b^2
52 cl_UDS_mul(b2_LSDptr,n2,a_LSDptr,n2,prod_LSDptr); // a*b^2
53 clear_loop_lsp(dest_LSDptr lspop n,n2-n);
54 shift1left_loop_lsp(dest_LSDptr,n+1); // (n+1 instead of n2 is ok)
55 subfrom_loop_lsp(prod_LSDptr,dest_LSDptr,n2);
60 // Bit complexity (N := len): O(M(N)).