src/modinteger/cl_MI_montgom.h

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