src/base/digitseq/cl_2DS_recip.cc

   1 // recip2adic().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_2DS.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_2D.h"
  13 #include "cl_DS.h"
  14
  15 // Break-even point of the Newton iteration vs. standard div2adic.
  16 const unsigned int recip2adic_threshold = 380;
  17
  18 void recip2adic (uintC len, const uintD* a_LSDptr, uintD* dest_LSDptr)
  19 {
  20         // Method:
  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).
  25         CL_ALLOCA_STACK;
  26         var uintL k = 0; // number of Newton steps
  27         var uintL n = len;
  28         while (n >= recip2adic_threshold) {
  29                 n = ceiling(n,2);
  30                 k++;
  31         }
  32         // Nonrecursive step.
  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);
  38         // Newton iteration.
  39         if (k > 0) {
  40                 var uintD* b2_LSDptr;
  41                 var uintD* prod_LSDptr;
  42                 num_stack_alloc(len+1,,b2_LSDptr=);
  43                 num_stack_alloc(2*(uintL)len,,prod_LSDptr=);
  44                 do {
  45                         // n = ceiling(len/2^k)
  46                         // Compute n2 = ceiling(len/2^(k-1)),
  47                         // then n = ceiling(n2/2).
  48                         k--;
  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);
  56                         n = n2;
  57                 } while (k > 0);
  58         }
  59 }
  60 // Bit complexity (N := len): O(M(N)).
  61