Montgomery Reduction

Purpose: Perform modular multiplication without any division, replacing % with shifts and multiplications. Typically 2-3× faster than % when the modulus is fixed and you do millions of multiplications — the standard speedup inside NTT, Miller-Rabin and Pollard’s rho.

The Idea

Fix an odd modulus $n$ and a radix $R=2^k>n$ (in code, $R=2^{32}$ or $2^{64}$, so division by $R$ is a free shift).

Work in the Montgomery domain: represent $a$ by $\bar{a}=aR\mathrm{mod}n$.

The reduction REDC(T) computes $T{R}^{-1}\mathrm{mod}n$ for $T<nR$:

$m=(T\mathrm{mod}R)\cdot n^{\prime }\mathrm{mod}R,t=\frac{T+mn}{R},$

where $n^{\prime }=-n^{-1}\mathrm{mod}R$. The numerator is divisible by $R$ by construction, so the division is exact and free. If $t\ge n$, subtract $n$ once.

Then $\overline{ab}=\mathrm{REDC}(\bar{a}\bar{b})$, so a domain multiplication is one REDC.

Code

struct Montgomery {
    uint64_t n, ninv, r2;                       // ninv = -n^{-1} mod 2^64
 
    Montgomery(uint64_t mod) : n(mod) {
        ninv = 1;
        for (int i = 0; i < 6; i++) ninv *= 2 - n * ninv;   // Newton, doubles bits
        ninv = -ninv;
        r2 = (__uint128_t(-(unsigned long long)n)) % n;      // R^2 mod n
    }
 
    uint64_t reduce(__uint128_t t) const {
        uint64_t m = (uint64_t)t * ninv;
        uint64_t res = (t + (__uint128_t)m * n) >> 64;
        return res >= n ? res - n : res;
    }
 
    uint64_t to(uint64_t a)   const { return reduce((__uint128_t)a * r2); }
    uint64_t from(uint64_t a) const { return reduce(a); }
    uint64_t mul(uint64_t a, uint64_t b) const { return reduce((__uint128_t)a * b); }
};

The Newton iteration ninv *= 2 - n * ninv doubles the number of correct bits each round, so six rounds suffice for 64 bits starting from 1 correct bit.

Complexity

• One REDC: two 64×64→128 multiplications plus a shift and a conditional subtract; no division
• Conversion in and out costs one REDC each — only worth it if you do many operations in the domain

Why It Works

$T+mn\equiv T(\mathrm{mod}n)$, and $m$ was chosen precisely so that $T+mn\equiv 0(\mathrm{mod}R)$. Hence $(T+mn)/R$ is an integer congruent to $T{R}^{-1}(\mathrm{mod}n)$. Bounding: $T<nR$ and $mn<Rn$, so $(T+mn)/R<2n$, meaning at most one conditional subtraction restores the range. ∎

Addition and subtraction are unchanged by the representation ($\overline{a+b}=\bar{a}+\bar{b}$), so only multiplication needs the reduction. Comparison and equality also work directly, since the map $a\mapsto aR\mathrm{mod}n$ is a bijection.

Odd modulus only

$n$ must be odd for $n^{-1}\mathrm{mod}{2}^k$ to exist. For even moduli, split with CRT or use Barrett reduction instead.

Variants / Use Cases

• Barrett reduction — works for even moduli and needs no domain conversion; slightly slower per operation but simpler to drop in
• NTT — the inner butterfly is a modmul; Montgomery is the standard optimisation
• Miller-Rabin and Pollard’s rho on 64-bit inputs
• RSA / ECC — Montgomery is why modular exponentiation is fast in every crypto library
• Modular Arithmetic — the topic page