Lucas' Theorem, in depth: binomial coefficients modulo a prime when n is larger than p, with proof, worked examples, code and the CRT extension

Why the usual factorial recipe breaks

The textbook formula is C(n, k) = n! / (k! (n - k)!). Modulo a prime p, division means multiplying by a modular inverse, and Fermat's little theorem gives the inverse of any a not divisible by p as a^(p-2) mod p. So the common implementation precomputes fact[i] = i! mod p and inv_fact[i] for every i up to n and returns fact[n] * inv_fact[k] * inv_fact[n-k] % p. With p = 1,000,000,007 and n up to a few million this is perfect, because no factorial up to n contains the factor p.

Now take p = 13 and n = 1000. Every factorial from 13! onwards is divisible by 13, so fact[1000] is 0 mod 13, and so are fact[300] and fact[700]. The formula becomes 0 times the inverse of 0, and 0 has no inverse. Fermat's trick happily returns 0^(11) = 0 as the "inverse", so the code prints 0 without any error, which is wrong: C(1000, 300) mod 13 is actually 10. The factors of p in the numerator and the denominator should cancel, but modular arithmetic has already thrown away the information needed to cancel them.

The theorem

Write n and k in base p: n = n_m p^m + ... + n_1 p + n_0 and k = k_m p^m + ... + k_1 p + k_0, where every digit is between 0 and p - 1 (pad the shorter number with leading zeros). Lucas' theorem, published in 1878, says that C(n, k) is congruent modulo p to the product of C(n_i, k_i) over all digit positions i.

Two things make this useful. First, every digit is smaller than p, so each small coefficient C(n_i, k_i) can be computed with the factorial recipe safely: no factorial below p contains p. Second, the convention C(a, b) = 0 when b is greater than a does real work. If any digit of k is larger than the matching digit of n, one factor is zero, so p divides C(n, k). The theorem turns one impossible computation into about log_p n tiny ones.

Why it is true

The proof rests on one fact: (1 + x)^p is congruent to 1 + x^p modulo p, as polynomials. Expanding (1 + x)^p gives coefficients C(p, j), and for 0 < j < p the coefficient C(p, j) = p! / (j! (p - j)!) has a p in the numerator and none in the denominator, so it is divisible by p. Only the first and last terms survive. Applying this repeatedly gives (1 + x)^(p^i) congruent to 1 + x^(p^i).

Now expand (1 + x)^n using the base-p digits of n. It equals the product over positions i of ((1 + x)^(p^i))^(n_i), which is congruent to the product of (1 + x^(p^i))^(n_i). Expanding each factor with the ordinary binomial theorem, the factor for position i contributes terms C(n_i, j_i) x^(j_i p^i) with 0 <= j_i <= n_i. Multiplying them, the coefficient of x^k collects the products of C(n_i, j_i) over all choices whose exponents sum to k. Because each j_i is at most n_i, which is below p, the sum of j_i p^i is a base-p representation, and base-p representations are unique. So exactly one choice reaches x^k, namely j_i = k_i, provided every k_i is at most n_i; otherwise no choice does and the coefficient is 0. Comparing with the coefficient of x^k on the left side, which is C(n, k), finishes the proof.

Worked example: C(1000, 300) modulo 13, 11 and 7

Start with p = 13. Since 13^2 = 169, divide: 1000 = 5 x 169 + 155 and 155 = 11 x 13 + 12, so 1000 has base-13 digits (5, 11, 12) from most significant to least. Similarly 300 = 1 x 169 + 131 and 131 = 10 x 13 + 1, giving (1, 10, 1). Pair the digits: C(5, 1) = 5, C(11, 10) = 11 and C(12, 1) = 12. Their product is 660, and 660 = 50 x 13 + 10, so C(1000, 300) mod 13 = 10, which matches comb(1000, 300) % 13.

Now p = 11. With 121 = 11^2, 1000 = 8 x 121 + 32 and 32 = 2 x 11 + 10, so the digits are (8, 2, 10). For 300, 300 = 2 x 121 + 58 and 58 = 5 x 11 + 3, so (2, 5, 3). The middle pair is C(2, 5), which is 0 because 5 exceeds 2. So 11 divides C(1000, 300), with no further arithmetic. Python agrees.

Finally p = 7. In base 7, 1000 is (2, 6, 2, 6) because 2 x 343 + 6 x 49 + 2 x 7 + 6 = 1000, and 300 is (0, 6, 0, 6). The factors are C(2, 0) = 1, C(6, 6) = 1, C(2, 0) = 1 and C(6, 6) = 1, so the answer is 1.

Implementation

The implementation has two parts: a factorial table of size p for small coefficients, and a loop that peels off one base-p digit of n and k per iteration. Building the inverse table backwards needs only one modular exponentiation.

from math import comb


class LucasBinomial:
    """C(n, k) mod p for a prime p and arbitrarily large n, k.

    Precomputation is O(p) time and memory, each query is O(log_p n).
    Practical while p is up to a few million.
    """

    def __init__(self, p: int):
        self.p = p
        self.fact = [1] * p
        for i in range(1, p):
            self.fact[i] = self.fact[i - 1] * i % p
        self.inv_fact = [1] * p
        self.inv_fact[p - 1] = pow(self.fact[p - 1], p - 2, p)  # Fermat: valid, fact[p-1] != 0
        for i in range(p - 1, 0, -1):
            self.inv_fact[i - 1] = self.inv_fact[i] * i % p

    def small(self, a: int, b: int) -> int:
        """C(a, b) mod p for 0 <= a, b < p."""
        if b > a:
            return 0
        return self.fact[a] * self.inv_fact[b] % self.p * self.inv_fact[a - b] % self.p

    def __call__(self, n: int, k: int) -> int:
        if k < 0 or k > n:
            return 0
        p, result = self.p, 1
        while n or k:
            a, b = n % p, k % p
            if b > a:
                return 0            # a borrow: p divides C(n, k)
            result = result * self.small(a, b) % p
            n //= p
            k //= p
        return result


if __name__ == "__main__":
    for p in (2, 3, 5, 7, 11, 13):
        lb = LucasBinomial(p)
        for n in range(0, 200):
            for k in range(0, n + 1):
                assert lb(n, k) == comb(n, k) % p, (p, n, k)
    print(LucasBinomial(13)(1000, 300))   # 10
    print(LucasBinomial(11)(1000, 300))   # 0
    print(LucasBinomial(7)(1000, 300))    # 1

The brute-force test at the bottom is not decoration. Off-by-one mistakes in the digit loop, a missing zero check or a table one element too short all pass a single hand-picked example and fail the exhaustive comparison immediately. In C++ or Java, use 64-bit integers and reduce after every multiplication.

Building the tables takes O(p) time and memory; each query takes O(log_p n) iterations, 17 for n = 10^18 and p = 13.

Parity and Kummer&#x27;s theorem

With p = 2 every digit is 0 or 1, and C(0, 1) = 0 is the only zero factor. So C(n, k) is odd exactly when every bit set in k is also set in n, which in code is (k & n) == k. The number of odd entries in row n of Pascal's triangle is therefore 2 raised to the number of set bits in n. Plot Pascal's triangle mod 2 and you get the Sierpinski triangle; Lucas' theorem is the reason.

Lucas tells you whether p divides C(n, k); Kummer's theorem (1852) tells you how many times. The exponent of p in C(n, k) equals the number of carries when k and n - k are added in base p. For p = 11 above, add 300 = (2, 5, 3) and 700 = (5, 8, 7) in base 11: 3 + 7 = 10 needs no carry, 5 + 8 = 13 carries, and 2 + 5 + 1 = 8 does not. One carry, so 11 divides C(1000, 300) exactly once. Kummer is what you need for questions like "what is the largest power of 2 dividing C(n, k)" and for prime-power moduli.

Composite moduli: CRT and prime powers

If the modulus m is a product of distinct primes, compute C(n, k) mod each prime with Lucas and combine the residues with the Chinese remainder theorem. For m = 1001 = 7 x 11 x 13 the residues from the worked example are 1, 0 and 10. CRT finds the unique x below 1001 with those residues, which is 946, and comb(1000, 300) % 1001 is indeed 946. The general tools, extended Euclid for the inverses and CRT itself, are covered in number theory for programmers.

def binom_mod_squarefree(n, k, primes):
    """C(n, k) mod (p1 * p2 * ...), primes distinct. Lucas per prime, then CRT."""
    m, x = 1, 0
    for p in primes:
        r = LucasBinomial(p)(n, k)
        # combine x (mod m) with r (mod p): find t with x + m*t = r (mod p)
        t = (r - x) * pow(m, -1, p) % p
        x, m = x + m * t, m * p
    return x


print(binom_mod_squarefree(1000, 300, [7, 11, 13]))   # 946, and comb(1000, 300) % 1001 == 946

If m contains a prime power such as 8 or 27, plain Lucas does not apply, because the digit factorization is only a congruence modulo p, not modulo p^e. The standard approach, published by Andrew Granville in 1997, computes n! with all factors of p removed, modulo p^e, using the periodicity of products of numbers coprime to p, and then restores the power of p that Kummer's theorem predicts. Factor m first with primes from a prime sieve, solve each prime power, then combine with CRT.

Choosing a method

Failure modes

Trade-offs

• Memory versus generality: tables are fastest but cost O(p) memory, which rules them out for large primes.
• Lucas versus exact arithmetic: for n of a few thousand, exact comb plus a modulus is simpler.

What to do next

1. Copy the LucasBinomial class and its exhaustive test, and run the test for primes 2 to 13.
2. Redo the C(1000, 300) example by hand for p = 5, then check it with the code.
3. Use Kummer's theorem to find the power of 2 dividing C(100, 50), then verify it with exact arithmetic.
4. Implement the CRT combination for a squarefree modulus of your choice, and read about Jacobi symbols in Legendre and Jacobi symbols to keep building your modular toolkit.
5. If your problem has prime-power moduli, read Granville's paper on binomial coefficients modulo prime powers before writing code.