The site's fourteenth Number Theory entry, and a genuinely different question from
every prior one: not a property of a single number, not a relationship between two, but a binomial
coefficient — C(n, k) mod p — where n and k can run to hundreds of
digits and p is a small, known prime. The obvious route computes n!,
k!, and (n−k)! mod p and combines them with
Fermat's-route modular inverse. That route works
right up until it doesn't: the moment k or n−k reaches p,
one of those factorials contains p itself as a factor, its value mod p
collapses to exactly 0, and 0 has no modular inverse — Fermat's route can't recover from that, and
(see Pitfalls) it doesn't fail loudly when it happens.
Lucas' Theorem sidesteps the whole problem: write n and k in base
p, then multiply the small binomial coefficient of each matching pair of digits —
C(ni, ki), always with ni, ki <
p — together mod p. No digit-level factorial ever reaches a multiple of
p, so the failure above never has a chance to happen. It's also the third entry to reuse
Modular Exponentiation's modPow
routine and the Fermat-inverse route directly,
extending that chain one link further — safely this time, because every digit fed into it is
guaranteed smaller than p.
Enter n, k, and a prime p, then compute
C(n, k) mod p digit by digit. The table shows each base-p digit of
n and k, the small binomial coefficient at that digit, and the running
product. Below it, three numbers are compared: Lucas' answer, a naive full-factorial-mod-p
route run on the whole n and k at once (skipped once the numbers are too
large to brute-force, which is itself the point), and an exact ground truth computed independently
with unbounded big-integer arithmetic (skipped at the same size). 7, 2, mod 5 is a
clean starting case. 10, 4, mod 3 is a real, correct zero. 7, 3, mod
3 is where the naive route breaks (k ≥ p). 6, 2, mod 4 uses
a composite modulus, which Lucas' Theorem requires be prime. 1015, 12345, mod
97 is far beyond where brute force could ever check — and still returns instantly.
Every other entry in this category reasons about one number's factorization or one modular
exponentiation. This one reasons about a whole generating function mod p. Start from a
fact provable directly from the binomial theorem: for a prime p, C(p, i) is
divisible by p for every 0 < i < p — its numerator p!
carries a factor of p that nothing in the denominator i!(p−i)! can
cancel, since every factor there is strictly smaller than p. That collapses
(1+x)p mod p down to just its first and last terms:
(1 + x)p ≡ 1 + xp (mod p)
— sometimes called the "freshman's dream," because it's the wrong shortcut every algebra student
tries with real numbers, and mod a prime it's actually true. Applying it i times gives
(1+x)pi ≡ 1 + xpi (mod p) for any
i. Now write n in base p as digits nm …
n0, so n = ∑ ni pi:
(1+x)n = ∏i (1+x)ni
pi ≡ ∏i (1 + xpi)ni
(mod p)
Expand each factor on the right with the ordinary binomial theorem: (1 +
xpi)ni = ∑j C(ni, j)
xj·pi, with 0 ≤ j ≤ ni < p. Picking
one term from each factor and multiplying them together produces x raised to
∑i ji·pi — and because every j is
itself smaller than p, that sum is literally the base-p representation of some
number, with digits ji. By uniqueness of base-p digit
representations, the only way that exponent can equal k is if ji =
ki at every digit position — no other combination of picks lands on x^k
at all. So the coefficient of xk on the right, which by construction is
C(n, k) mod p on the left, is exactly ∏i C(ni,
ki) mod p — Lucas' Theorem, read straight off both sides of the same equation.
function modPow(base, exp, mod) { // identical to Modular Exponentiation's own routine
base = ((base % mod) + mod) % mod;
let result = 1n % mod;
while (exp > 0n) {
if (exp & 1n) result = (result * base) % mod;
base = (base * base) % mod;
exp >>= 1n;
}
return result;
}
function modInverseFermat(a, p) { // identical to the Fermat-inverse route; p MUST be prime
return modPow(a, p - 2n, p);
}
// n, k both strictly less than p here — the whole point of decomposing into digits first
function smallBinomial(n, k, p) {
if (k < 0n || k > n) return 0n;
let num = 1n, den = 1n;
for (let i = 0n; i < k; i++) {
num = (num * ((n - i) % p)) % p;
den = (den * ((i + 1n) % p)) % p;
}
return (num * modInverseFermat(den, p)) % p;
}
function lucas(n, k, p) { // p MUST be prime
if (k < 0n) return 0n;
let result = 1n, nn = n, kk = k;
while (nn > 0n || kk > 0n) { // BUG (see Pitfalls): must check BOTH nn and kk
const ni = nn % p, ki = kk % p;
if (ki > ni) return 0n; // BUG (see Pitfalls): omitting this is a real, checked bug elsewhere in this family
result = (result * smallBinomial(ni, ki, p)) % p;
nn /= p;
kk /= p;
}
return result;
}Two routines this site already built, called on numbers that never exceed p — that's
the entire algorithm. The loop terminates once both n's and k's digits are
exhausted, multiplying one small binomial coefficient in per iteration.
The naive route this page opens with doesn't fail loudly when k or n−k reaches p — it
silently returns 0, right only by coincidence. Run the ordinary multiplicative
formula for C(n, k) mod p directly on the whole n and k (no
digit decomposition at all), and the running denominator product 1 × 2 × …
× k is guaranteed to include a multiple of p the moment k ≥
p — at that point the denominator is congruent to 0 mod p, and
modInverseFermat(0, p) (which is modPow(0, p−2, p)) returns exactly 0,
not an error. The formula then multiplies by that 0 and returns 0 for every such input,
regardless of the true answer. A worked case: n = 7, k = 3, p = 3 (k ≥
p) — the naive route returns 0, but C(7,3) = 35 and
35 mod 3 = 2, the value Lucas' Theorem gives. Forcing k ≥ p across
30,000 random trials, the naive route returned exactly 0 in 100.0% of them (checked
directly) — and against the true value, that 0 was wrong 27.2% of the time
(8,152 of 30,000). The other 72.8% happen to have a genuinely zero answer too (a real, if
coincidental, agreement — C(n,k) mod p = 0 is common whenever adding k and
n−k in base p needs a carry), which is exactly what makes this failure
mode dangerous: it looks indistinguishable from a correct 0 nearly three times out of four.
A composite modulus breaks the per-digit modular inverse the same way it breaks Fermat's
route on its own page — silently, with a plausible-looking wrong number. Lucas' Theorem
requires p prime; nothing in the code above checks that. Worked case:
n = 6, k = 2, p = 4 (composite). Base-4 digits: n = 124,
k = 024. The high digit is fine (C(1,0) = 1), but the low digit
needs C(2, 2) mod 4 — and computing that via
modInverseFermat(2, 4) asks for the inverse of 2 mod 4, which doesn't exist
(gcd(2, 4) = 2 ≠ 1). modPow(2, 4−2, 4) = modPow(2, 2, 4) = 0 comes
back anyway — not an error, just a wrong number — so that digit's term is computed as 0
instead of the true C(2,2) = 1, and the whole product comes out 0. The
true value is C(6,2) mod 4 = 15 mod 4 = 3. Swept across 30,000 random trials over
composite moduli 4–16, Lucas' Theorem's own routine disagreed with the true value in
10.9% of them (3,256 of 30,000) — lower than the composite-modulus failure rate
measured on the plain Fermat-inverse page, because many digit pairs here happen to avoid hitting a
non-invertible denominator at all, not because the underlying problem is any less real.
The loop must check both n's and k's remaining digits, not just n's — dropping the second
check silently misses trailing digits of k whenever k > n. When k > n, the
true answer is always 0 (there's no way to choose more items than exist). Lucas' Theorem
still gets there correctly, but only because while (nn > 0n || kk > 0n) keeps
looping until both numbers run out of digits — a version that loops only
while (nn > 0n) stops as soon as n's digits are exhausted, even if
k still has nonzero digits above that point that were never checked against
ki > ni = 0. Worked case: n = 5, k = 100, p =
3. The correct routine returns 0 (as it must, since 100 >
5); the shortened-loop version returns 2 — a wrong, nonzero-looking
answer. Forcing k > n (so the true answer is always 0) across 30,000 random trials,
the shortened-loop version returned a nonzero, wrong value in 8.2% of them (2,468 of
30,000); the rest still landed on 0 by coincidence, whenever a mismatched digit happened to fall
within the digits the shortened loop still checks.
O(p logp n) modular multiplications: each digit's
smallBinomial call does up to O(p) work (its multiplicative formula runs
ki < p times), and there are O(logp n) digits to
process. Each digit also pays one O(log p) call into modInverseFermat, which
this bound absorbs since p dominates log p for essentially every prime this
algorithm would actually be used with. Checked directly: 100,000 random trials against an
independent exact big-integer ground truth (unbounded BigInt arithmetic, no modulus
until the final reduction), 0 mismatches.
The cost depends on p, not on the size of C(n, k) itself, which is the
entire point — lucas(1015, 12345n, 97n) and lucas(10n, 3n, 97n)
run at essentially the same speed, while the exact value of C(1015, 12345)
itself has tens of thousands of digits and the naive full-factorial-mod-p route above would need to
process a denominator product 12,345 terms long before ever reaching the digit decomposition this
algorithm uses instead.
This site's guide, Choosing a Number Theory Algorithm, compares this entry against the other thirteen Number Theory entries side by side.