Chinese Remainder Theorem, in depth: constructive proof, Garner mixed radix, non-coprime merging, multi-prime NTT and RSA-CRT

The theorem and why it holds

Let m1, m2, ..., mk be pairwise coprime positive integers and M their product. For any residues r1, ..., rk the system x = ri (mod mi) for every i has exactly one solution with 0 <= x < M. Equivalently, the map from x mod M to the tuple (x mod m1, ..., x mod mk) is a bijection that respects addition and multiplication. In algebra this is the ring isomorphism Z/MZ to the product of the Z/miZ.

Uniqueness is short. If x and y both solve the system, every mi divides x - y; since the moduli are coprime, their product M divides x - y, so x and y agree modulo M. Existence is constructive. Let Mi = M / mi. Because Mi shares no factor with mi, it has an inverse yi modulo mi, which the extended Euclidean algorithm finds. Then x = sum of ri Mi yi, reduced mod M, works: for term j not equal to i, Mj is divisible by mi and vanishes, and term i gives ri times 1.

The isomorphism view is what makes CRT useful in engineering. Arithmetic on one huge modulus becomes independent arithmetic on several small ones, each fitting a machine word, with no carries between them.

A worked example by hand

Take the classic puzzle: a number leaves remainder 2 when divided by 3, remainder 3 when divided by 5 and remainder 2 when divided by 7. Here M = 105. Using the textbook formula: M1 = 35, and 35 mod 3 = 2 whose inverse mod 3 is 2. M2 = 21, 21 mod 5 = 1, inverse 1. M3 = 15, 15 mod 7 = 1, inverse 1. So x = 2 x 35 x 2 + 3 x 21 x 1 + 2 x 15 x 1 = 140 + 63 + 30 = 233, and 233 mod 105 = 23. Check: 23 = 7 x 3 + 2, 23 = 4 x 5 + 3, 23 = 3 x 7 + 2.

Notice the intermediate value 233 exceeds M. With three small moduli nobody cares, but with many 60-bit moduli the textbook formula multiplies M-sized numbers k times. That is the motivation for Garner.

Garner mixed-radix reconstruction

Garner's algorithm computes x in mixed-radix form: x = v0 + v1 m1 + v2 m1 m2 + ..., where each digit vi lies in [0, mi). Each digit is found modulo mi alone, so every step is a word-size subtraction and multiplication. Only the final assembly, if you need it, touches big integers, and often you do not: if the goal is x modulo some other number p, evaluate the mixed-radix sum mod p and never build x at all.

from math import gcd, prod

def garner(residues, moduli):
    """Mixed-radix digits v with x = v0 + v1*m0 + v2*m0*m1 + ...
    Every intermediate value stays below the largest modulus."""
    k = len(moduli)
    # inv[j][i] = m_j^{-1} mod m_i, precomputable once per modulus set
    inv = [[pow(moduli[j], -1, moduli[i]) if j < i else 0 for i in range(k)]
           for j in range(k)]
    v = []
    for i in range(k):
        t = residues[i] % moduli[i]
        for j in range(i):
            t = (t - v[j]) * inv[j][i] % moduli[i]
        v.append(t)
    return v

def from_mixed_radix(v, moduli):
    x, base = 0, 1
    for digit, m in zip(v, moduli):
        x += digit * base
        base *= m
    return x

def crt_coprime(residues, moduli):
    return from_mixed_radix(garner(residues, moduli), moduli)

assert garner([2, 3, 2], [3, 5, 7]) == [2, 2, 1]
assert crt_coprime([2, 3, 2], [3, 5, 7]) == 23

Trace the example. v0 = 2. For m = 5: (3 - 2) x inverse of 3 mod 5, which is 2, gives v1 = 2. For m = 7: (2 - 2) x inverse of 3 mod 7 = 0, then (0 - 2) x inverse of 5 mod 7, which is 3, gives -6 = 1 mod 7. So x = 2 + 2 x 3 + 1 x 15 = 23. The inverse table depends only on the moduli, so libraries compute it once and reuse it for millions of reconstructions; the per-call work is then k(k-1)/2 multiply-reduce steps.

When the moduli are not coprime

Real inputs do not always have coprime moduli. Scheduling problems, cycle detection and many contest tasks give congruences like x = 2 (mod 6) and x = 8 (mod 10). The rule: a solution exists if and only if every pair of residues agrees modulo the gcd of their moduli, and when it exists it is unique modulo the least common multiple. Here gcd is 2, and 2 and 8 are both even, so a solution exists modulo 30.

Merge two congruences at a time. Write x = r1 + m1 t and require m1 t = r2 - r1 (mod m2). Divide the whole congruence by g = gcd(m1, m2), which is legal exactly when g divides r2 - r1, and invert m1/g modulo m2/g.

def crt_merge(r1, m1, r2, m2):
    """Combine x = r1 (mod m1) and x = r2 (mod m2) for ANY moduli.
    Returns (r, lcm) or None if the congruences contradict."""
    g = gcd(m1, m2)
    if (r2 - r1) % g:
        return None
    m2g = m2 // g
    # solve m1 * t = r2 - r1 (mod m2) by dividing through by g
    t = ((r2 - r1) // g * pow(m1 // g, -1, m2g)) % m2g if m2g > 1 else 0
    lcm = m1 // g * m2
    return (r1 + m1 * t) % lcm, lcm

def crt_general(residues, moduli):
    r, m = 0, 1
    for ri, mi in zip(residues, moduli):
        merged = crt_merge(r, m, ri, mi)
        if merged is None:
            return None
        r, m = merged
    return r, m

assert crt_merge(2, 6, 8, 10) == (8, 30)
assert crt_merge(1, 6, 2, 4) is None
assert crt_general([2, 3, 2], [3, 5, 7]) == (23, 105)

For x = 1 (mod 6) and x = 2 (mod 4), the gcd is 2 but 2 - 1 is odd: the first says x is odd and the second says x is even, so no solution exists. Returning None rather than a wrong number is the most important line in the function. Pairwise merging also handles the coprime case, at the cost of big-integer intermediates, so it is the safe default when inputs are untrusted.

Multi-prime NTT and residue number systems

Fast polynomial multiplication with the number-theoretic transform works modulo a prime of the form c x 2^k + 1, such as 998244353 = 119 x 2^23 + 1. If you need the exact integer convolution of two arrays with coefficients up to 10^9 and length 10^5, each output coefficient can reach about 10^23, far above any single prime. The fix is to run the convolution modulo three NTT-friendly primes and reconstruct each coefficient with CRT, since the product of the three primes, about 7.9 x 10^25, exceeds the largest possible value.

P = (998244353, 167772161, 469762049)  # c * 2^k + 1 primes, all with root 3
# conv_mod(a, b, p) is any NTT convolution modulo prime p
def big_convolution(a, b):
    lanes = [conv_mod(a, b, p) for p in P]
    out = []
    for r0, r1, r2 in zip(*lanes):
        out.append(crt_coprime([r0, r1, r2], list(P)))
    return out

The correctness condition is easy to forget: the product of the moduli must exceed the true maximum value, including sign if coefficients can be negative, in which case map results above M/2 back to negatives. The same lane idea appears in the floating-point FFT world as splitting coefficients into high and low halves; CRT achieves exactness without rounding analysis.

Residue number systems push the idea into hardware and cryptography. Lattice-based homomorphic encryption schemes such as BFV and CKKS represent a huge ciphertext modulus as a product of machine-word primes and keep every coefficient in residue form, so additions and multiplications run lane by lane on CPUs and GPUs. The hard operations in RNS are the ones that need the whole number at once, such as comparison, division and rescaling, which require either partial reconstruction or base-conversion tricks. That is the general trade: CRT makes ring operations embarrassingly parallel and makes ordering operations expensive.

RSA-CRT and the fault attack

RSA private operations are the best-known production use. Instead of computing m^d mod N with N = pq, a library computes m^dp mod p and m^dq mod q with exponents reduced modulo p - 1 and q - 1 by Fermat, then joins the two halves with a two-modulus Garner step. Modular exponentiation costs roughly the cube of the modulus length, so two half-size exponentiations cost about a quarter of one full one; real libraries report roughly three to four times faster decryption and signing. That is why private keys store dp, dq and the inverse of q mod p alongside d, as described in the RSA key format.

def rsa_crt_sign(m, p, q, d, e):
    dp, dq = d % (p - 1), d % (q - 1)
    q_inv = pow(q, -1, p)
    s1 = pow(m, dp, p)             # half-size exponent, half-size modulus
    s2 = pow(m, dq, q)
    h = (q_inv * (s1 - s2)) % p    # Garner with two moduli
    s = s2 + h * q
    if pow(s, e, p * q) != m % (p * q):   # verify before release
        raise RuntimeError("fault detected; signature withheld")
    return s

The final check is not decoration. In 1997 Boneh, DeMillo and Lipton showed that if a fault, such as a voltage glitch, a cosmic-ray bit flip or a buggy multiplier, corrupts only one half, the released signature s' is correct modulo q and wrong modulo p. Then s'^e - m is divisible by q but not p, and gcd(s'^e - m, N) reveals q, and with it the whole key, from a single bad signature. Verifying with the cheap public exponent before release, or comparing against a recomputation, is the standard countermeasure, and mature libraries do it.

Failure modes

Most CRT bugs fall into a few groups. Assuming coprimality when it does not hold makes pow raise an error in Python and returns garbage in hand-written inverse code in other languages; validate with gcd or use the general merge. Overflow in fixed-width languages: (t - v[j]) * inv can exceed 64 bits when moduli are near 2^63, so use 128-bit intermediates or Montgomery multiplication. Negative residues: in C and Java the % operator can return negatives, so normalise. Insufficient modulus product in multi-prime convolution silently wraps results. And in cryptography, skipping fault verification or using non-constant-time arithmetic on secret residues turns a speed-up into a key leak.

Trade-offs

Textbook reconstruction is short and fine for a handful of small moduli. Garner needs a precomputed inverse table but keeps intermediates word-sized and can stop early to produce x mod p directly. Pairwise merging handles arbitrary moduli and detects contradictions at the cost of growing intermediates. In systems design, RNS buys carry-free parallel arithmetic and pays for it whenever you need comparison, division or overflow detection, so it suits workloads dominated by multiply-accumulate, which is exactly what polynomial and lattice arithmetic are.

What to do next

1. Solve the 3, 5, 7 example by hand with both the textbook formula and Garner, and confirm 23.
2. Copy the Python above and add a randomised test comparing garner against brute force for small moduli.
3. Extend crt_merge tests with contradictory pairs and moduli where one divides the other.
4. Implement a three-prime NTT convolution and check it against Python big-integer multiplication.
5. Read your RSA library source and find where it verifies CRT signatures before release.
6. If you write C or Rust, rewrite garner with 128-bit intermediates and fuzz it near 2^63.