ii. Barrett and Montgomery Reduction

% by a run-time divisor costs 20–90 cycles; the compiler only turns division by a compile-time constant into multiplications. When a hot loop reduces modulo a value read from input (or a 64-bit modulus in Miller–Rabin / Pollard rho), replace the division by multiplications:

Tool	Modulus	Idea
Barrett	any $m$, run time	precompute $\lfloor 2^{64}/m\rfloor$; the estimated quotient is off by at most 1
Lemire fastmod	$m<2^{32}$, operands $<2^{32}$	64-bit fractional reciprocal, two multiplications
Montgomery	odd $m<2^{62}$	keep $aR\mathrm{mod}m$; reduce products with REDC (no division at all)
mulmod64	$m\le 7.2\cdot 10^{18}$	long double estimate of $\lfloor ab/m\rfloor$, exact correction

using u64 = unsigned long long; using u128 = __uint128_t;
struct Barrett {                                           // x % m for any x < 2^64
    u64 m, im;
    Barrett(u64 m) : m(m), im(~0ULL / m) {}
    u64 reduce(u64 x) const { u64 q = (u64)(((u128)im * x) >> 64), r = x - q * m; return r >= m ? r - m : r; }
    u64 mul(u64 a, u64 b) const { return reduce(a * b); } // a, b < m < 2^32
};
struct FastMod32 {                                         // Lemire: x % d for x, d < 2^32
    u64 M; unsigned d;
    FastMod32(unsigned d) : M(~0ULL / d + 1), d(d) {}
    unsigned reduce(unsigned x) const { u64 low = M * x; return (unsigned)(((u128)low * d) >> 64); }
};
struct Montgomery64 {                                      // odd n < 2^62, values kept as a * 2^64 mod n
    u64 n, ninv, r2;
    Montgomery64(u64 n) : n(n), ninv(n) {
        for (int i = 0; i < 6; i++) ninv *= 2 - n * ninv;  // n^{-1} mod 2^64 (Newton: bits double)
        ninv = -ninv; r2 = (u64)(-(u128)n % n);           // -n^{-1}, and 2^128 mod n
    }
    u64 redc(u128 t) const { u64 m = (u64)t * ninv; u64 r = (u64)((t + (u128)m * n) >> 64); return r >= n ? r - n : r; }
    u64 to(u64 a) const { return redc((u128)(a % n) * r2); }
    u64 val(u64 a) const { return redc(a); }
    u64 mul(u64 a, u64 b) const { return redc((u128)a * b); }   // both in Montgomery form
    u64 add(u64 a, u64 b) const { a += b; return a >= n ? a - n : a; }
    u64 pow(u64 a, u64 e) const { u64 r = to(1); for (; e; e >>= 1, a = mul(a, a)) if (e & 1) r = mul(r, a); return r; }
};
u64 mulmod64(u64 a, u64 b, u64 m) {                        // a, b < m <= 7.2e18, no __int128
    ll r = (ll)(a * b - m * (u64)(1.L / m * a * b));
    return r + m * (r < 0) - m * (r >= (ll)m);
}

Montgomery form is closed under addition and subtraction, equality is preserved (the map is a bijection), and only multiplication needs redc. Convert in once, compute many products (Miller–Rabin powers, Pollard rho iterations, NTT butterflies, RSA/ECC exponentiation — Montgomery is why modular exponentiation is fast in every crypto library), convert out once. In practice Barrett is roughly 4x faster than % by a run-time divisor, and Montgomery another 2-3x faster than Barrett when the modulus is fixed for many multiplications — the gap is this large because Barrett still does two full-width multiplications per %, while Montgomery’s redc replaces the division entirely and the conversion overhead amortizes away.

Montgomery needs an odd modulus

$n$ must be odd for $n^{-1}\mathrm{mod}{2}^k$ to exist (the Newton iteration requires it). For an even modulus, use Barrett reduction instead, or split the even part out via CRT (CRT) and run Montgomery on the odd factor.

Proof: Barrett and Lemire

Barrett: let $\text{im}=\lfloor (2^{64}-1)/m\rfloor$, so $\frac{2^{64}}{m}-1-\frac{1}{m}<\text{im}\le \frac{2^{64}}{m}$. Then $q=\lfloor \text{im}\cdot x/2^{64}\rfloor \le x/m$, and $\frac{x}{m}-q<1+\frac{x(1+1/m)}{2^{64}}<2$ for $x<2^{64}$. So $q\in \{\lfloor x/m\rfloor -1,\lfloor x/m\rfloor \}$ and $r=x-qm\in [0,2m)$: one conditional subtraction gives $x\mathrm{mod}m$. (mul needs $ab<2^{64}$.)

Lemire: $M=\lceil 2^{64}/d\rceil$ (for $d$ not a power of two; exact otherwise), and $\text{low}=(Mx)\mathrm{mod}{2}^{64}$ is the fractional part of $x/d$ scaled by $2^{64}$, up to an error $<x\cdot \frac{Md-2^{64}}{d}\cdot \frac{1}{1}\le x<2^{32}$ units of $2^{-64}$. Multiplying the fractional part by $d$ and taking the integer part recovers $x\mathrm{mod}d$ exactly because the accumulated error $x(Md-2^{64})<2^{32}\cdot 2^{32}=2^{64}$ never reaches the next integer of the scaled product.

Proof: Montgomery REDC and the Newton inverse

REDC: with $R=2^{64}$ and $n^{\prime }=-n^{-1}\mathrm{mod}R$, set $m=(t\mathrm{mod}R)n^{\prime }\mathrm{mod}R$. Then $t+mn\equiv t-t\equiv 0(\mathrm{mod}R)$, so $u=(t+mn)/R$ is an integer, and $uR\equiv t(\mathrm{mod}n)$ gives $u\equiv t{R}^{-1}(\mathrm{mod}n)$. If $t<nR$ then $u<(nR+Rn)/R=2n$, so one subtraction suffices. With inputs $\bar{a}=aR$, $\bar{b}=bR$ (both $<n$), $\text{redc}(\bar{a}\bar{b})=abR=\overline{ab}$; to multiplies by $R^2$ and reduces once; val reduces once to remove the factor $R$. Newton: if $nx\equiv 1(\mathrm{mod}{2}^k)$, then $x^{\prime }=x(2-nx)$ satisfies $1-n{x}^{\prime }=(1-nx{)}^2\equiv 0(\mathrm{mod}{2}^{2k})$; starting from $x=n$ ($n^2\equiv 1(\mathrm{mod}8)$, 3 correct bits) six iterations reach 64 bits. Range: $n<2^{62}$ keeps $t+mn<2^{128}$.

Proof: the long double mulmod (for $m<2^{62}$)

Let $c=\lfloor ab/m\rfloor <m$. Each of the three long double operations (64-bit mantissa) has relative error at most $2^{-64}$, so the computed value of $ab/m$ is within about $3\cdot 2^{-64}\cdot 2^{62}<1$ of the truth, and after truncation $q\in \{c-1,c,c+1\}$. Hence $ab-qm\in [-m,2m)$. The expression a * b - m * q is evaluated modulo $2^{64}$ but its true value lies in $[-m,2m)\subset [-2^{63},2^{63})$, so the cast to signed gives it exactly, and one correction in either direction lands in $[0,m)$. KACTL’s finer error analysis extends the range to $m\le 7.2\cdot 10^{18}$ (verified by random testing); near $2^{63}$ it fails, so use __int128 there. It requires the x87 80-bit long double.