src/modinteger/cl_MI_montgom.h

   1 // m > 1 odd, Montgomery representation
   2
   3 // We use Montgomery's modular multiplication trick
   4 // [Peter L. Montgomery: Modular multiplication without trial division,
   5 //  Mathematics of Computation 44 (1985), 519-521.]
   6 //
   7 // Assume we want to compute modulo M, M odd. V and N will be chosen
   8 // so that V*N==1 mod M and that (a,b) --> a*b*V mod M can be more easily
   9 // computed than (a,b) --> a*b mod M. Then, we have a ring isomorphism
  10 //   (Z/MZ, +, * mod M)  \isomorph  (Z/MZ, +, (a,b) --> a*b*V mod M)
  11 //   x mod M             -------->  x*N mod M
  12 // It is thus preferrable to use x*N mod M as a "representation" of x mod M,
  13 // especially for computations which involve at least several multiplications.
  14 //
  15 // The precise algorithm to compute a*b*V mod M, given a and b, and the choice
  16 // of N and V depend on M and on the hardware. The general idea is this:
  17 // Choose N = 2^n, so that division by N is easy. Recall that V == N^-1 mod M.
  18 // 1. Given a and b as m-bit numbers (M <= 2^m), compute a*b in full
  19 //    precision.
  20 // 2. Write a*b = c*N+d, i.e. split it into components c and d.
  21 // 3. Now a*b*V = c*N*V+d*V == c+d*V mod M.
  22 // 4. Instead of computing d*V mod M
  23 //    a. by full multiplication and then division mod M, or
  24 //    b. by left shifts: repeated application of
  25 //          x := 2*x+(0 or 1); if (x >= M) { x := x-M; }
  26 //    we compute
  27 //    c. by right shifts (recall that d*V == d*2^-n mod M): repeated application
  28 //       of   if (x odd) { x := (x+M)/2; } else { x := x/2; }
  29 // Usually one will choose N = 2^m, so that c and d have both m bits.
  30 // Several variations are possible: In step 4 one can implement the right
  31 // shifts in hardware. Or (for example when N = 2^160 and working on a
  32 // 32-bit machine) one can do 32 shift steps at the same time:
  33 // Choose M' == M^-1 mod 2^32 and compute n/32 times
  34 //       x := (x - ((x mod 2^32) * M' mod 2^32) * M) / 2^32.
  35
  36 // Here, we choose to use Montgomery representation only if |V| can be chosen
  37 // to be very small, and in that case we compute d*V mod M using standard
  38 // multiplication and division.
  39 // So we choose N = 2^n with 0 < n <= m (the larger n, the better) and hope
  40 // that it will yield V with |V| < 2^k. We thus replace a division of
  41 // 2m bits by m bits (cost: approx. m^2) by a multiplication of n bits with
  42 // k bits (cost: approx. n*k) and a division of max(2m-n,n+k) bits by m bits
  43 // (cost: approx. (max(2m-n,n+k)-m)*m). Of course, U*M+V*N=1 implies (roughly)
  44 // n+k >= m. It is worth it when
  45 //         m^2 > n*k + (n+k-m)*m   and   m^2 > n*k + (m-n)*m
  46 // <==>  3*m^2 > (n+m)*(k+m)       and     m > k
  47 // <==   3/2*m > k+m                (assume n to be near m)
  48 // <==>    m/2 > k .
  49 //
  50 // How to find N and V:
  51 // U*M+V*N=1 means that U = (M mod 2^n)^-1 = U_m mod 2^n, where
  52 // U_m := (M mod 2^m)^-1 (2-adic reciprocal). |V| < 2^(m/2) is more or less
  53 // equivalent to |V*N| < 2^(n+m/2) <==> |U|*M < 2^(n+m/2) <==> |U| < n-m/2
  54 // <==> the most significant m/2 bits of |U| are all equal. So we search
  55 // for a bit string of at least m/2+1 equal bits in U_m, which has m bits.
  56 // Very easy: take the middle bit of U_m, look how many bits adjacent to it
  57 // (at the left and at the right) have the same value. Say these are the
  58 // bits n-1,...,n-l. (Choose n and l as large as possible. k = m-l + O(1).)
  59 // If l < m/2, forget it. Else fix n and compute V = (1-U*M)/2^n.
  60 //
  61 // It is now clear that only very few moduli M will allow such a good
  62 // choice of N and V, but in these cases the Montgomery multiplication
  63 // reduces the multiplication complexity by a large constant factor.
  64
  65
  66 class cl_heap_modint_ring_montgom : public cl_heap_modint_ring {
  67         SUBCLASS_cl_heap_modint_ring()
  68 public:
  69         // Constructor.
  70         cl_heap_modint_ring_montgom (const cl_I& M, uintL m, uintL n, const cl_I& V);
  71         // Virtual destructor.
  72         ~cl_heap_modint_ring_montgom () {}
  73         // Additional information.
  74         uintL m; // M = 2^m
  75         uintL n; // N = 2^n, n <= m
  76         cl_I V;
  77 };
  78
  79 // Assuming 0 <= x < 2^(2m), return  V*x mod M.
  80 static inline const cl_I montgom_redc (cl_heap_modint_ring_montgom* R, const cl_I& x)
  81 {
  82         return mod((x >> R->n) + (R->V * ldb(x,cl_byte(R->n,0))), R->modulus);
  83 }
  84
  85 static const _cl_MI montgom_canonhom (cl_heap_modint_ring* _R, const cl_I& x)
  86 {
  87         var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
  88         return _cl_MI(R, mod(x << R->n, R->modulus));
  89 }
  90
  91 static const cl_I montgom_retract (cl_heap_modint_ring* _R, const _cl_MI& x)
  92 {
  93         var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
  94         return montgom_redc(R,x.rep);
  95 }
  96
  97 static const _cl_MI montgom_one (cl_heap_modint_ring* _R)
  98 {
  99         var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
 100         var cl_I zr = (cl_I)1 << R->n;
 101         return _cl_MI(R, R->n == R->m ? zr - R->modulus : zr);
 102 }
 103
 104 static const _cl_MI montgom_mul (cl_heap_modint_ring* _R, const _cl_MI& x, const _cl_MI& y)
 105 {
 106         var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
 107         return _cl_MI(R, montgom_redc(R,x.rep * y.rep));
 108 }
 109
 110 static const _cl_MI montgom_square (cl_heap_modint_ring* _R, const _cl_MI& x)
 111 {
 112         var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
 113         return _cl_MI(R, montgom_redc(R,square(x.rep)));
 114 }
 115
 116 static const cl_MI_x montgom_recip (cl_heap_modint_ring* _R, const _cl_MI& x)
 117 {
 118         var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
 119         var const cl_I& xr = x.rep;
 120         var cl_I u, v;
 121         var cl_I g = xgcd(xr,R->modulus,&u,&v);
 122         // g = gcd(x,M) = x*u+M*v
 123         if (eq(g,1))
 124                 return cl_MI(R, mod((minusp(u) ? u + R->modulus : u) << (2*R->n), R->modulus));
 125         if (zerop(xr))
 126                 cl_error_division_by_0();
 127         return cl_notify_composite(R,xr);
 128 }
 129
 130 static const cl_MI_x montgom_div (cl_heap_modint_ring* _R, const _cl_MI& x, const _cl_MI& y)
 131 {
 132         var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
 133         var const cl_I& yr = y.rep;
 134         var cl_I u, v;
 135         var cl_I g = xgcd(yr,R->modulus,&u,&v);
 136         // g = gcd(y,M) = y*u+M*v
 137         if (eq(g,1))
 138                 return cl_MI(R, mod((x.rep * (minusp(u) ? u + R->modulus : u)) << R->n, R->modulus));
 139         if (zerop(yr))
 140                 cl_error_division_by_0();
 141         return cl_notify_composite(R,yr);
 142 }
 143
 144 #define montgom_addops std_addops
 145 static cl_modint_mulops montgom_mulops = {
 146         montgom_one,
 147         montgom_canonhom,
 148         montgom_mul,
 149         montgom_square,
 150         std_expt_pos,
 151         montgom_recip,
 152         montgom_div,
 153         std_expt,
 154         std_reduce_modulo,
 155         montgom_retract
 156 };
 157
 158 // Constructor.
 159 inline cl_heap_modint_ring_montgom::cl_heap_modint_ring_montgom (const cl_I& M, uintL _m, uintL _n, const cl_I& _V)
 160         : cl_heap_modint_ring (M, &std_setops, &montgom_addops, &montgom_mulops),
 161           m (_m), n (_n), V (_V)
 162 {}
 163
 164 static cl_heap_modint_ring* try_make_modint_ring_montgom (const cl_I& M)
 165 {
 166         if (!oddp(M))
 167                 return NULL;
 168         var uintL m = integer_length(M);
 169         CL_ALLOCA_STACK;
 170         var uintC len;
 171         var const uintD* M_LSDptr;
 172         I_to_NDS_nocopy(M, ,len=,M_LSDptr=,cl_false,);
 173         if (lspref(M_LSDptr,len-1)==0) { len--; } // normalize
 174         // Compute U as 2-adic inverse of M.
 175         var uintD* U_LSDptr;
 176         num_stack_alloc(len,,U_LSDptr=);
 177         recip2adic(len,M_LSDptr,U_LSDptr);
 178         // Look at U's bits.
 179         #define U_bit(i) (lspref(U_LSDptr,floor(i,intDsize)) & ((uintD)1 << ((i)%intDsize)))
 180         var uintL i_min;
 181         var uintL i_max;
 182         var uintL i = floor(m,2);
 183         var cl_boolean negative;
 184         if (U_bit(i)) {
 185                 for (; --i > 0; )
 186                         if (!U_bit(i)) break;
 187                 i_min = i+1;
 188                 i = floor(m,2);
 189                 for (; ++i < m; )
 190                         if (!U_bit(i)) break;
 191                 i_max = i;
 192                 negative = cl_true;
 193         } else {
 194                 for (; --i > 0; )
 195                         if (U_bit(i)) break;
 196                 i_min = i+1;
 197                 i = floor(m,2);
 198                 for (; ++i < m; )
 199                         if (U_bit(i)) break;
 200                 i_max = i;
 201                 negative = cl_false;
 202         }
 203         #undef U_bit
 204         // OK, all the bits i_max-1..i_min of U are equal.
 205         if (i_max - i_min <= floor(m,2))
 206                 return NULL;
 207         var uintL n = i_max;
 208         // Turn U (mod 2^n) into a signed integer.
 209         if (n % intDsize) {
 210                 if (negative)
 211                         lspref(U_LSDptr,floor(n,intDsize)) |= (uintD)(-1) << (n % intDsize);
 212                 else
 213                         lspref(U_LSDptr,floor(n,intDsize)) &= ((uintD)1 << (n % intDsize)) - 1;
 214         }
 215         var uintL U_len = ceiling(n,intDsize);
 216         var cl_I U = DS_to_I(U_LSDptr lspop U_len,U_len);
 217         var cl_I V_N = 1 - U*M;
 218         if (ldb_test(V_N,cl_byte(n,0)))
 219                 cl_abort();
 220         var cl_I V = V_N >> n;
 221         return new cl_heap_modint_ring_montgom(M,m,n,V);
 222 }