Combinatorics, in depth: counting rules, binomials modulo a prime, Lucas's theorem, inclusion-exclusion, Catalan numbers and ranking combinations

The sum and product rules

Everything starts with two rules. The sum rule: if a choice can be made in one of several disjoint ways, add the counts. The product rule: if a choice is a sequence of independent steps, multiply the counts per step. A grid search over 4 learning rates, 3 batch sizes, 5 warmup schedules and 2 optimisers has 4 x 3 x 5 x 2 = 120 configurations by the product rule.

Most counting mistakes come from applying these rules to cases that are not disjoint or not independent. The rest of combinatorics is largely machinery for correcting those two situations: division to remove orderings you did not mean to count, and inclusion-exclusion to remove overlaps.

The four selection cases

Choosing k items from n distinct items splits on two questions: does order matter, and may an item be chosen more than once?

The unordered, no-repetition case divides the ordered count by k!, because each subset appears in k! orders. C(52, 5) = 2,598,960 poker hands. The unordered-with-repetition case is stars and bars: distributing k identical items into n bins is arranging k stars and n-1 bars in a row, so the count is C(n+k-1, n-1). Splitting 10 identical GPU-hours across 3 jobs gives C(12, 2) = 66 ways.

A worked example that combines rules: a feature-selection sweep picks exactly 2 of 6 candidate features and, independently, one of 3 regularisation strengths for each of 2 model heads. The product rule gives C(6, 2) x 3² = 15 x 9 = 135 runs. If each run takes 20 minutes, that is 45 GPU-hours, which you now know before you launch.

Exact binomials and overflow

Python integers are unbounded, so math.comb(n, k) is exact. In languages with fixed-width integers, binomials overflow quickly. C(66, 33) = 7,219,428,434,016,265,740 is the largest central binomial that fits in a signed 64-bit integer; C(67, 33) needs an unsigned one, and C(68, 34) overflows both. The multiplicative formula avoids factorials but must be ordered so every intermediate division is exact:

def binom(n, k):
    if k < 0 or k > n:
        return 0
    k = min(k, n - k)
    r = 1
    for i in range(k):
        r = r * (n - i) // (i + 1)   # exact: r is C(n, i+1) after this step
    return r

After step i the running value is C(n, i+1), an integer, so the floor division never truncates. In C or Java, the product r * (n - i) can still overflow before the division; divide by a gcd first or use 128-bit intermediates. For many queries with small n, Pascal's rule C(n, k) = C(n-1, k-1) + C(n-1, k) fills a table with additions only, at O(n²) memory, as in Subset Sum counting tables.

Binomials modulo a prime

Contest problems and many probabilistic computations ask for counts modulo a large prime such as 10⁹+7. Division is not defined directly in modular arithmetic, but when p is prime every nonzero value a has an inverse a^p-2 by Fermat's little theorem. Precompute factorials and inverse factorials once; then each binomial is three multiplications.

MOD = 10**9 + 7

def build(n, mod=MOD):
    fact = [1] * (n + 1)
    for i in range(1, n + 1):
        fact[i] = fact[i - 1] * i % mod
    inv = [1] * (n + 1)
    inv[n] = pow(fact[n], mod - 2, mod)      # one modular exponentiation
    for i in range(n, 0, -1):
        inv[i - 1] = inv[i] * i % mod        # 1/(i-1)! = i / i!
    return fact, inv

fact, inv = build(10**6)

def nCr(n, r):
    if r < 0 or r > n:
        return 0
    return fact[n] * inv[r] % MOD * inv[n - r] % MOD

print(nCr(10, 3))      # 120
print(nCr(1000, 500))  # 159835829, matches math.comb(1000, 500) % MOD

Building the table is O(n) with a single exponentiation, a trick explained in Modular Exponentiation, in depth. The method has a hard precondition: n must be smaller than p, because otherwise fact[n] contains the factor p and is zero.

When n reaches p: Lucas's theorem

When the prime is small or n is huge, use Lucas's theorem: write n and k in base p, and C(n, k) mod p is the product of C(n_i, k_i) mod p over the digits, with any digit pair where k_i > n_i making the whole result zero.

def lucas(n, k, p):
    f = [1] * p
    for i in range(1, p):
        f[i] = f[i - 1] * i % p
    def small(a, b):
        if b > a:
            return 0
        return f[a] * pow(f[b], p - 2, p) * pow(f[a - b], p - 2, p) % p
    res = 1
    while n or k:
        res = res * small(n % p, k % p) % p
        n //= p
        k //= p
    return res

Example: C(20, 6) mod 7. In base 7, 20 is digits (2, 6) and 6 is (0, 6), so the result is C(2, 0) x C(6, 6) = 1. The true value is 38,760, and 38,760 mod 7 is indeed 1. The factorial-table method returns 0 here because 20! is divisible by 7, a silent wrong answer rather than an error. For a composite modulus, factor it, compute modulo each prime power, and combine with the Chinese remainder theorem from Number Theory Algorithms; prime powers need an extension of Lucas that tracks the factors of p separately.

Inclusion-exclusion: surjections and derangements

Inclusion-exclusion counts objects that avoid every one of several bad events A₁..A_m: start from the total, subtract objects in each single event, add back those in each pair, subtract triples, and so on. It works whenever intersections are easy to count even though the union is not.

Surjections. How many ways can 5 distinct jobs be assigned to 3 workers so every worker gets at least one? Let A_j be assignments where worker j gets nothing. The answer is the sum over j of (-1)^j C(3, j) (3-j)⁵ = 243 - 96 + 3 - 0 = 150. Dividing by 3! when workers are interchangeable gives the Stirling number of the second kind S(5, 3) = 25, the number of ways to partition 5 items into 3 non-empty groups.

Derangements. Permutations with no fixed point, such as a secret-santa draw where nobody gets their own name, count as D(n) = n! times the sum of (-1)ⁱ/i!. A cleaner recurrence is D(n) = (n-1)(D(n-1) + D(n-2)), giving 0, 1, 2, 9, 44, 265, 1854 for n = 1 to 7. D(n)/n! tends to 1/e, about 0.368, which is why a random shuffle leaves at least one item in place about 63 percent of the time.

from math import comb

def surjections(n, k):
    return sum((-1)**j * comb(k, j) * (k - j)**n for j in range(k + 1))

def derangements(n):
    if n == 0:
        return 1
    a, b = 1, 0                      # D(0), D(1)
    for i in range(2, n + 1):
        a, b = b, (i - 1) * (a + b)
    return b

print(surjections(5, 3), derangements(5))   # 150 44

Catalan and Stirling numbers

Some sequences recur so often they deserve names. The Catalan numbers C_n = C(2n, n)/(n+1), which run 1, 1, 2, 5, 14, 42, 132, 429, count balanced strings of n bracket pairs, binary tree shapes with n nodes, ways to triangulate a polygon with n+2 sides, and monotone grid paths that never cross the diagonal. The common structure is a first-return decomposition: split at the point where the first bracket closes, giving the recurrence C_n+1 = sum of C_i C_n-i. Whenever a recursive decomposition splits a structure into an independent left and right part, check the first few counts against this sequence before deriving anything.

Stirling numbers of the second kind satisfy S(n, k) = k S(n-1, k) + S(n-1, k-1): the last item either joins one of k existing groups or starts its own. That recurrence fills a table in O(nk) and is the dynamic-programming view of the same count inclusion-exclusion gave above; see Dynamic programming, in depth for the general technique.

Ranking and unranking combinations

Sometimes you need the i-th combination directly: to shard a search space across workers, to resume an enumeration after a crash, or to sample a uniform random subset by drawing an index. The combinatorial number system maps each k-subset {c₁ < ... < c_k} of {0..n-1} to the integer sum of C(c_i, i), a bijection onto 0..C(n, k)-1.

from math import comb

def rank(c):                     # c sorted ascending
    return sum(comb(x, i + 1) for i, x in enumerate(c))

def unrank(N, k):
    out = []
    for i in range(k, 0, -1):
        x = i - 1
        while comb(x + 1, i) <= N:   # largest x with comb(x, i) <= N
            x += 1
        out.append(x)
        N -= comb(x, i)
    return out[::-1]

print(rank([1, 3, 4]), unrank(8, 3))   # 8 [1, 3, 4]

For 3-subsets of 5 items the ranks run 0 to 9, and [1, 3, 4] has rank C(1,1) + C(3,2) + C(4,3) = 8. Giving worker w the ranks from w x C(n,k)/W up to the next worker's start divides the space evenly without any coordination. The linear scan in unrank is fine for small n; replace it with binary search when n is large.

Failure modes

Trade-offs

Exact big-integer arithmetic is simplest and correct, but its cost grows with the size of the numbers; C(10⁶, 5 x 10⁵) has about 300,000 decimal digits. Modular counts are fast and fixed-size but lose magnitude, so you cannot compare them or turn them into probabilities. Log-space estimates via lgamma are ideal for deciding whether an enumeration is feasible and useless for exact answers. Pascal tables trade O(n²) memory for division-free arithmetic under any modulus, including composites. When counts must be convolved, for example the number of ways to reach each total across several independent choices, polynomial multiplication via the Fast Fourier Transform turns an O(n²) combination step into O(n log n).

What to do next

1. For each counting problem, write down whether order matters and whether repetition is allowed before choosing a formula.
2. Verify every formula on a tiny case by brute-force enumeration with itertools before trusting it on large inputs.
3. Implement the factorial and inverse-factorial table and check it against math.comb modulo the prime.
4. Add a guard that rejects n at least p, and keep a Lucas implementation ready for small moduli.
5. Practise inclusion-exclusion on surjections and derangements until the sign pattern is automatic.
6. Use ranking and unranking the next time you need to split a combination space across workers.