Counting: nCr modulo p and Catalan numbers

C(n, k), read “n choose k”, counts the ways to choose k items from n when the order does not matter: C(n, k) = n! / (k! × (n − k)!). Its values grow fast, and C(1000, 500) already has 300 digits, so programs usually compute it modulo a prime p such as 10⁹ + 7. Modulo a prime, build tables of n! and of the inverses of the factorials once, in O(n), and every query is then two multiplications. When the modulus is not prime, Pascal’s rule C(n, k) = C(n − 1, k − 1) + C(n − 1, k) uses additions only and works for any modulus, in O(n × k). The Catalan numbers, C(2n, n) / (n + 1) = 1, 1, 2, 5, 14, 42 and so on, count the shapes of binary search trees, balanced strings of brackets and many other structures.

What C(n, k) counts

Choosing k items is the same as leaving out the other n − k, so C(n, k) = C(n, n − k), and there is no way to choose more items than there are, so C(n, k) = 0 when k > n:

import math
from itertools import combinations

print("3 of the 5 letters a to e, order ignored:")
print(" ", " ".join("".join(c) for c in combinations("abcde", 3)))
print(f"  that is C(5, 3) = {math.comb(5, 3)}, and C(5, 2) = {math.comb(5, 2)}: keeping 3 is leaving out 2")
print(f"6 numbers of 49: C(49, 6) = {math.comb(49, 6):,}")
print(f"C(1000, 500) has {len(str(math.comb(1000, 500)))} digits")
print(f"C(4, 7) = {math.comb(4, 7)}: there is no way to choose 7 items from 4")

3 of the 5 letters a to e, order ignored:
  abc abd abe acd ace ade bcd bce bde cde
  that is C(5, 3) = 10, and C(5, 2) = 10: keeping 3 is leaving out 2
6 numbers of 49: C(49, 6) = 13,983,816
C(1000, 500) has 300 digits
C(4, 7) = 0: there is no way to choose 7 items from 4

Python’s math.comb returns the exact value at any size, which makes it the reference to test modular code against. Exact values are expensive, though: C(1000, 500) has 300 digits, and a program that answers a hundred thousand such queries works with remainders instead.

Factorial tables modulo a prime

Division modulo p means multiplying by an inverse, so C(n, k) mod p = n! × (k!)⁻¹ × ((n − k)!)⁻¹ mod p. Building the tables takes two passes. The first goes up: fact[i] = fact[i − 1] × i mod p. The second goes down from a single inverse: inv_fact[N] = fact[N]^(p − 2) by Fermat’s little theorem, and then inv_fact[i − 1] = inv_fact[i] × i, since 1 / (i − 1)! = i / i!. One fast exponentiation and 2N multiplications build both tables.

import math

P = 10**9 + 7
N = 1_000_000

# fact[i] = i! mod P, built upwards; inv_fact[i] = (i!)^-1 mod P, built downwards from one Fermat inverse.
fact = [1] * (N + 1)
for i in range(1, N + 1):
    fact[i] = fact[i - 1] * i % P
inv_fact = [1] * (N + 1)
inv_fact[N] = pow(fact[N], P - 2, P)
for i in range(N, 0, -1):
    inv_fact[i - 1] = inv_fact[i] * i % P


def comb_mod(n, k):
    if k < 0 or k > n:
        return 0
    return fact[n] * inv_fact[k] % P * inv_fact[n - k] % P


def comb_by_product(n, k):
    """An independent check: n(n-1)...(n-k+1) / k!, with a single inverse at the end."""
    top = bottom = 1
    for i in range(k):
        top = top * (n - i) % P
        bottom = bottom * (i + 1) % P
    return top * pow(bottom, P - 2, P) % P


print(f"C(49, 6) mod P = {comb_mod(49, 6):,}")
print(f"C(100000, 50000) mod P = {comb_mod(10**5, 5 * 10**4)}, exact math.comb agrees: {comb_mod(10**5, 5 * 10**4) == math.comb(10**5, 5 * 10**4) % P}")
print(f"C(1000000, 500000) mod P = {comb_mod(N, N // 2)}, the product formula agrees: {comb_mod(N, N // 2) == comb_by_product(N, N // 2)}")
print(f"C(10, 11) = {comb_mod(10, 11)}, C(10, -1) = {comb_mod(10, -1)}")
print(f"each query: 2 multiplications and 2 reductions after building {2 * N:,} table entries")

C(49, 6) mod P = 13,983,816
C(100000, 50000) mod P = 149033233, exact math.comb agrees: True
C(1000000, 500000) mod P = 996692777, the product formula agrees: True
C(10, 11) = 0, C(10, -1) = 0
each query: 2 multiplications and 2 reductions after building 2,000,000 table entries

Each query is now two multiplications, C(49, 6) comes back as 13,983,816, and the values past the modulus agree with exact arithmetic for C(100,000, 50,000) and with a separate product formula for C(1,000,000, 500,000). Queries with k < 0 or k > n must return 0 before any table lookup.

Pascal’s triangle for any modulus

Every entry of Pascal’s triangle is the sum of the two above it, because a set of k items from n either contains the last item, leaving k − 1 to choose from the first n − 1, or does not, leaving all k to choose from them:

That rule needs no division, so it works modulo any number. Factorial tables fail modulo 10⁹ = 2⁹ × 5⁹, because every factorial from 40! on is a multiple of 10⁹ and has no inverse:

import math

M = 10**9  # not a prime: 10^9 = 2^9 × 5^9

print(f"40! mod 10^9 = {math.factorial(40) % M}: factorials are multiples of 10^9 from 40! on, so no inverse exists")

row = [1]  # row n of Pascal's triangle, modulo M: C(n, k) = C(n - 1, k - 1) + C(n - 1, k)
additions = 0
for n in range(1, 1001):
    row = [1] + [(row[k - 1] + row[k]) % M for k in range(1, n)] + [1]
    additions += n - 1
print(f"row 1000 built with {additions:,} additions")
print(f"C(1000, 500) mod 10^9 = {row[500]}")
print(f"every entry of row 1000 matches math.comb(1000, k) % 10^9: {all(row[k] == math.comb(1000, k) % M for k in range(1001))}")

40! mod 10^9 = 0: factorials are multiples of 10^9 from 40! on, so no inverse exists
row 1000 built with 499,500 additions
C(1000, 500) mod 10^9 = 821216320
every entry of row 1000 matches math.comb(1000, k) % 10^9: True

Row 1000 took 499,500 additions and matches math.comb(1000, k) % 10**9 in every entry. A table of rows up to n is O(n²) work and memory, fine for n of a few thousand.

Catalan numbers

Count the binary search trees that can hold the keys 1 to n. Choosing key r as the root leaves r − 1 keys for the left subtree and n − r for the right one, and any left shape goes with any right shape, so the counts satisfy Cₙ = Σ C(r − 1) × C(n − r) over every root r, starting from C₀ = 1. These are the Catalan numbers, and their closed form is C(2n, n) / (n + 1):

import math
from functools import lru_cache


def shapes_by_dp(n):
    """count[i] = the number of BST shapes with i keys: pick a root, split the other keys left and right."""
    count = [1] + [0] * n
    for i in range(1, n + 1):
        count[i] = sum(count[left] * count[i - 1 - left] for left in range(i))
    return count


@lru_cache(maxsize=None)
def all_shapes(n):
    """Every BST shape with n keys, as nested tuples (left, right), for a brute-force count."""
    if n == 0:
        return [None]
    return [(l, r) for left in range(n) for l in all_shapes(left) for r in all_shapes(n - 1 - left)]


dp = shapes_by_dp(10)
formula = [math.comb(2 * n, n) // (n + 1) for n in range(11)]
print(" n  shapes by DP  C(2n, n)/(n + 1)  listed one by one")
for n in range(1, 11):
    listed = f"{len(all_shapes(n)):,}" if n <= 8 else "-"
    print(f"{n:>2}  {dp[n]:>12,}  {formula[n]:>16,}  {listed:>17}")

 n  shapes by DP  C(2n, n)/(n + 1)  listed one by one
 1             1                 1                  1
 2             2                 2                  2
 3             5                 5                  5
 4            14                14                 14
 5            42                42                 42
 6           132               132                132
 7           429               429                429
 8         1,430             1,430              1,430
 9         4,862             4,862                  -
10        16,796            16,796                  -

The DP, the formula and a brute-force listing agree: 14 trees for 4 keys and 16,796 for 10. The same numbers count the balanced strings of n pairs of brackets, the ways to cut a convex polygon with n + 2 corners into triangles, and the grid paths that never cross a diagonal.

Lucas’s theorem computes C(n, k) modulo a small prime p for huge n by splitting n and k into base-p digits.

Exercise and self-check

Answer many C(n, k) queries modulo 10⁹ + 7, building the tables once.

def ncr_mod(queries):
    """C(n, k) mod 10**9 + 7 for each (n, k): 0 when k < 0 or k > n."""
    raise NotImplementedError("write ncr_mod")

import math
import random

from ncr import ncr_mod

P = 10**9 + 7


def test_small_values():
    """small binomial coefficients, below the modulus"""
    assert ncr_mod([(5, 2), (49, 6), (10, 0), (10, 10), (0, 0)]) == [10, 13_983_816, 1, 1, 1]


def test_k_outside_the_range():
    """k below 0 or above n gives 0"""
    assert ncr_mod([(4, 7), (4, -1), (0, 1), (3, 4)]) == [0, 0, 0, 0]


def test_symmetry():
    """C(n, k) = C(n, n - k)"""
    pairs = [(1000, 3), (1000, 997), (777, 300), (777, 477)]
    a, b, c, d = ncr_mod(pairs)
    assert a == b and c == d


def test_values_past_the_modulus():
    """C(1000, 500) and C(10^6, 5 × 10^5), whose exact values are far larger than P"""
    assert ncr_mod([(1000, 500)]) == [math.comb(1000, 500) % P]
    assert ncr_mod([(10**6, 5 * 10**5)]) == [996_692_777]


def test_ten_thousand_random_queries():
    """10,000 random queries with n up to 3,000, against exact arithmetic"""
    rng = random.Random(17)
    queries = []
    for _ in range(10_000):
        n = rng.randint(0, 3000)
        queries.append((n, rng.randint(-5, n + 5)))
    answers = ncr_mod(queries)
    assert sum(answers) == 4_836_647_905_907
    for (n, k), answer in list(zip(queries, answers))[:200]:
        assert answer == (math.comb(n, k) % P if 0 <= k <= n else 0), (n, k)

/** C(n, k) mod 1,000,000,007 for each [n, k]: 0 when k < 0 or k > n. */
export function ncrMod(queries) {
  throw new Error('write ncrMod');
}

import { test, assert } from 'mysmartcopilot:test';
import { ncrMod } from './ncr.mjs';

/** A small seeded generator (Park and Miller), so every run builds the same queries. */
function generator(seed) {
  const random = () => (seed = (seed * 48271) % 2147483647) / 2147483647;
  return (k) => Math.floor(random() * k);
}

test('small binomial coefficients, below the modulus', () => {
  assert.deepEqual(ncrMod([[5, 2], [49, 6], [10, 0], [10, 10], [0, 0]]), [10, 13983816, 1, 1, 1]);
});

test('k below 0 or above n gives 0', () => {
  assert.deepEqual(ncrMod([[4, 7], [4, -1], [0, 1], [3, 4]]), [0, 0, 0, 0]);
});

test('C(n, k) = C(n, n - k)', () => {
  const [a, b, c, d] = ncrMod([[1000, 3], [1000, 997], [777, 300], [777, 477]]);
  assert.ok(a === b && c === d);
});

test('C(1000, 500) and C(10^6, 5 × 10^5), whose exact values are far larger than P', () => {
  assert.deepEqual(ncrMod([[1000, 500]]), [159835829]);
  assert.deepEqual(ncrMod([[1000000, 500000]]), [996692777]);
});

test('10,000 random queries with n up to 3,000: their answers add up to a known total', () => {
  const pick = generator(17);
  const queries = [];
  for (let i = 0; i < 10000; i++) {
    const n = pick(3001);
    queries.push([n, pick(n + 11) - 5]);
  }
  const answers = ncrMod(queries);
  assert.equal(answers.length, 10000);
  assert.equal(answers.reduce((s, x) => s + x, 0), 4837030980937);
});

Why the tables are right, and nCr in five languages Pro: opens with a Pro, Premium or Ultimate pass

Catalan numbers proved two ways, and found in four problems Pro: opens with a Pro, Premium or Ultimate pass

Lucas's theorem when n is at least p, and a practice set Pro: opens with a Pro, Premium or Ultimate pass

Interview questions

Key takeaways

• C(n, k) = n! / (k! (n − k)!) = C(n, n − k), and it is 0 when k < 0 or k > n.
• Modulo a prime p > n, build fact and inv_fact once, with a single Fermat inverse, and answer each query with two multiplications.
• Pascal’s rule uses additions only, so it works for any modulus, at O(n²) for all rows up to n.
• Catalan numbers, C(2n, n) / (n + 1), count BST shapes, balanced brackets, triangulations and paths that stay below a diagonal.