Montgomery Multiplication

Modular multiplications x y mod n show up constantly in number theory and crypto (primality, factoring, RSA). A naive product-then-divide is slow: division by n dominates, especially for big integers.

Montgomery multiplication replaces that division by adding a multiple of n that clears the low bits, then discarding those bits (cheap when the auxiliary modulus r is a power of two).

This lesson follows CP-Algorithms — Montgomery multiplication.

Montgomery representation

Fix odd n (typical) and r ≥ n with gcd(n,r)=1. In practice r = 2ᵐ so mod// by r are shifts and masks.

x̄ := x · r  (mod n)

Entering the space is itself a modular multiply (expensive once). After that, add / subtract / equality / gcd-with-n work as usual, because they respect the common factor r.

Multiplication does not:

x̄ · ȳ ≡ (x y) · r²  (mod n)

but we want

x̄ * ȳ := (x y) · r  (mod n)

So define

x̄ * ȳ := x̄ · ȳ · r⁻¹  (mod n)

Montgomery reduction (REDC)

Need x · r⁻¹ mod n efficiently. Since gcd(n,r)=1, Bézout gives r⁻¹, n' with

r · r⁻¹ + n · n' = 1

(n' ≡ n⁻¹ (mod r)). Classic REDC:

function reduce(x):          # assume 0 ≤ x < n·r
    q = (x mod r) · n' mod r
    a = (x − q·n) / r
    if a < 0: a += n
    return a

No division by n — only ops mod// r. Final correction is at most one add of n.

Applications:

Montgomery mul: reduce(x̄ · ȳ)
Leave the space: reduce(x̄) recovers x

Fast inverse mod `2ᵏ` (Newton)

a·x ≡ 1 (mod 2ᵏ)  ⇒  a·x·(2 − a·x) ≡ 1 (mod 2²ᵏ)

Start from x=1 (inverse mod 2) and double correct bits each squaring-style step. For r=2¹²⁸, about 7 iterations suffice.

Implementation notes (128-bit)

GCC provides __int128 / __uint128, so 64×64 → mod n is easy. For 128×128, split into high/low 64-bit limbs (u256::mult), then REDC with r=2¹²⁸.

struct Montgomery {
    Montgomery(u128 n) : mod(n), inv(1) {
        for (int i = 0; i < 7; i++)
            inv *= 2 - n * inv;   // wraps mod 2^128
    }
    u128 reduce(u256 x);          // q = x.low * inv; use high half
    u128 mult(u128 a, u128 b) { return reduce(u256::mult(a, b)); }
    u128 mod, inv;
};

Fast transform into the space

x̄ = x · r ≡ x * r²   (Montgomery mul by r²)

Precompute r2 = r² mod n, then init(x) = mult(x, r2). One construction: start from r ≡ −n (mod n), build 2⁴ in Montgomery form, square five times to get 2¹²⁸ ≡ r in the space (= r² in ordinary residues).

🧪 Toy REDC (n=97, r=128)

Check: Montgomery

Quiz