* All Files have been modified for inclusion of namespace cln;

[cln.git] / src / modinteger / cl_MI_montgom.h

// m > 1 odd, Montgomery representation

namespace cln {

// We use Montgomery's modular multiplication trick
// [Peter L. Montgomery: Modular multiplication without trial division,
//  Mathematics of Computation 44 (1985), 519-521.]
//
// Assume we want to compute modulo M, M odd. V and N will be chosen
// so that V*N==1 mod M and that (a,b) --> a*b*V mod M can be more easily
// computed than (a,b) --> a*b mod M. Then, we have a ring isomorphism
//   (Z/MZ, +, * mod M)  \isomorph  (Z/MZ, +, (a,b) --> a*b*V mod M)
//   x mod M             -------->  x*N mod M
// It is thus preferrable to use x*N mod M as a "representation" of x mod M,
// especially for computations which involve at least several multiplications.
//
// The precise algorithm to compute a*b*V mod M, given a and b, and the choice
// of N and V depend on M and on the hardware. The general idea is this:
// Choose N = 2^n, so that division by N is easy. Recall that V == N^-1 mod M.
// 1. Given a and b as m-bit numbers (M <= 2^m), compute a*b in full
//    precision.
// 2. Write a*b = c*N+d, i.e. split it into components c and d.
// 3. Now a*b*V = c*N*V+d*V == c+d*V mod M.
// 4. Instead of computing d*V mod M
//    a. by full multiplication and then division mod M, or
//    b. by left shifts: repeated application of
//          x := 2*x+(0 or 1); if (x >= M) { x := x-M; }
//    we compute
//    c. by right shifts (recall that d*V == d*2^-n mod M): repeated application
//       of   if (x odd) { x := (x+M)/2; } else { x := x/2; }
// Usually one will choose N = 2^m, so that c and d have both m bits.
// Several variations are possible: In step 4 one can implement the right
// shifts in hardware. Or (for example when N = 2^160 and working on a
// 32-bit machine) one can do 32 shift steps at the same time:
// Choose M' == M^-1 mod 2^32 and compute n/32 times
//       x := (x - ((x mod 2^32) * M' mod 2^32) * M) / 2^32.

// Here, we choose to use Montgomery representation only if |V| can be chosen
// to be very small, and in that case we compute d*V mod M using standard
// multiplication and division.
// So we choose N = 2^n with 0 < n <= m (the larger n, the better) and hope
// that it will yield V with |V| < 2^k. We thus replace a division of
// 2m bits by m bits (cost: approx. m^2) by a multiplication of n bits with
// k bits (cost: approx. n*k) and a division of max(2m-n,n+k) bits by m bits
// (cost: approx. (max(2m-n,n+k)-m)*m). Of course, U*M+V*N=1 implies (roughly)
// n+k >= m. It is worth it when
//         m^2 > n*k + (n+k-m)*m   and   m^2 > n*k + (m-n)*m
// <==>  3*m^2 > (n+m)*(k+m)       and     m > k
// <==   3/2*m > k+m                (assume n to be near m)
// <==>    m/2 > k .
//
// How to find N and V:
// U*M+V*N=1 means that U = (M mod 2^n)^-1 = U_m mod 2^n, where
// U_m := (M mod 2^m)^-1 (2-adic reciprocal). |V| < 2^(m/2) is more or less
// equivalent to |V*N| < 2^(n+m/2) <==> |U|*M < 2^(n+m/2) <==> |U| < n-m/2
// <==> the most significant m/2 bits of |U| are all equal. So we search
// for a bit string of at least m/2+1 equal bits in U_m, which has m bits.
// Very easy: take the middle bit of U_m, look how many bits adjacent to it
// (at the left and at the right) have the same value. Say these are the
// bits n-1,...,n-l. (Choose n and l as large as possible. k = m-l + O(1).)
// If l < m/2, forget it. Else fix n and compute V = (1-U*M)/2^n.
//
// It is now clear that only very few moduli M will allow such a good
// choice of N and V, but in these cases the Montgomery multiplication
// reduces the multiplication complexity by a large constant factor.


class cl_heap_modint_ring_montgom : public cl_heap_modint_ring {
        SUBCLASS_cl_heap_modint_ring()
public:
        // Constructor.
        cl_heap_modint_ring_montgom (const cl_I& M, uintL m, uintL n, const cl_I& V);
        // Virtual destructor.
        ~cl_heap_modint_ring_montgom () {}
        // Additional information.
        uintL m; // M = 2^m
        uintL n; // N = 2^n, n <= m
        cl_I V;
};

// Assuming 0 <= x < 2^(2m), return  V*x mod M.
static inline const cl_I montgom_redc (cl_heap_modint_ring_montgom* R, const cl_I& x)
{
        return mod((x >> R->n) + (R->V * ldb(x,cl_byte(R->n,0))), R->modulus);
}

static const _cl_MI montgom_canonhom (cl_heap_modint_ring* _R, const cl_I& x)
{
        var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
        return _cl_MI(R, mod(x << R->n, R->modulus));
}

static const cl_I montgom_retract (cl_heap_modint_ring* _R, const _cl_MI& x)
{
        var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
        return montgom_redc(R,x.rep);
}

static const _cl_MI montgom_one (cl_heap_modint_ring* _R)
{
        var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
        var cl_I zr = (cl_I)1 << R->n;
        return _cl_MI(R, R->n == R->m ? zr - R->modulus : zr);
}

static const _cl_MI montgom_mul (cl_heap_modint_ring* _R, const _cl_MI& x, const _cl_MI& y)
{
        var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
        return _cl_MI(R, montgom_redc(R,x.rep * y.rep));
}

static const _cl_MI montgom_square (cl_heap_modint_ring* _R, const _cl_MI& x)
{
        var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
        return _cl_MI(R, montgom_redc(R,square(x.rep)));
}

static const cl_MI_x montgom_recip (cl_heap_modint_ring* _R, const _cl_MI& x)
{
        var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
        var const cl_I& xr = x.rep;
        var cl_I u, v;
        var cl_I g = xgcd(xr,R->modulus,&u,&v);
        // g = gcd(x,M) = x*u+M*v
        if (eq(g,1))
                return cl_MI(R, mod((minusp(u) ? u + R->modulus : u) << (2*R->n), R->modulus));
        if (zerop(xr))
                cl_error_division_by_0();
        return cl_notify_composite(R,xr);
}

static const cl_MI_x montgom_div (cl_heap_modint_ring* _R, const _cl_MI& x, const _cl_MI& y)
{
        var cl_heap_modint_ring_montgom* R = (cl_heap_modint_ring_montgom*)_R;
        var const cl_I& yr = y.rep;
        var cl_I u, v;
        var cl_I g = xgcd(yr,R->modulus,&u,&v);
        // g = gcd(y,M) = y*u+M*v
        if (eq(g,1))
                return cl_MI(R, mod((x.rep * (minusp(u) ? u + R->modulus : u)) << R->n, R->modulus));
        if (zerop(yr))
                cl_error_division_by_0();
        return cl_notify_composite(R,yr);
}

#define montgom_addops std_addops
static cl_modint_mulops montgom_mulops = {
        montgom_one,
        montgom_canonhom,
        montgom_mul,
        montgom_square,
        std_expt_pos,
        montgom_recip,
        montgom_div,
        std_expt,
        std_reduce_modulo,
        montgom_retract
};

// Constructor.
inline cl_heap_modint_ring_montgom::cl_heap_modint_ring_montgom (const cl_I& M, uintL _m, uintL _n, const cl_I& _V)
        : cl_heap_modint_ring (M, &std_setops, &montgom_addops, &montgom_mulops),
          m (_m), n (_n), V (_V)
{}

static cl_heap_modint_ring* try_make_modint_ring_montgom (const cl_I& M)
{
        if (!oddp(M))
                return NULL;
        var uintL m = integer_length(M);
        CL_ALLOCA_STACK;
        var uintC len;
        var const uintD* M_LSDptr;
        I_to_NDS_nocopy(M, ,len=,M_LSDptr=,cl_false,);
        if (lspref(M_LSDptr,len-1)==0) { len--; } // normalize
        // Compute U as 2-adic inverse of M.
        var uintD* U_LSDptr;
        num_stack_alloc(len,,U_LSDptr=);
        recip2adic(len,M_LSDptr,U_LSDptr);
        // Look at U's bits.
        #define U_bit(i) (lspref(U_LSDptr,floor(i,intDsize)) & ((uintD)1 << ((i)%intDsize)))
        var uintL i_min;
        var uintL i_max;
        var uintL i = floor(m,2);
        var cl_boolean negative;
        if (U_bit(i)) {
                for (; --i > 0; )
                        if (!U_bit(i)) break;
                i_min = i+1;
                i = floor(m,2);
                for (; ++i < m; )
                        if (!U_bit(i)) break;
                i_max = i;
                negative = cl_true;
        } else {
                for (; --i > 0; )
                        if (U_bit(i)) break;
                i_min = i+1;
                i = floor(m,2);
                for (; ++i < m; )
                        if (U_bit(i)) break;
                i_max = i;
                negative = cl_false;
        }
        #undef U_bit
        // OK, all the bits i_max-1..i_min of U are equal.
        if (i_max - i_min <= floor(m,2))
                return NULL;
        var uintL n = i_max;
        // Turn U (mod 2^n) into a signed integer.
        if (n % intDsize) {
                if (negative)
                        lspref(U_LSDptr,floor(n,intDsize)) |= (uintD)(-1) << (n % intDsize);
                else
                        lspref(U_LSDptr,floor(n,intDsize)) &= ((uintD)1 << (n % intDsize)) - 1;
        }
        var uintL U_len = ceiling(n,intDsize);
        var cl_I U = DS_to_I(U_LSDptr lspop U_len,U_len);
        var cl_I V_N = 1 - U*M;
        if (ldb_test(V_N,cl_byte(n,0)))
                cl_abort();
        var cl_I V = V_N >> n;
        return new cl_heap_modint_ring_montgom(M,m,n,V);
}

}  // namespace cln