Binomial Coefficients: Precomputation & Modular Arithmetic

June 2026 · 20 min read

Knowing the formula C(n, r) = n! / (r!(n−r)!) is one thing; computing it fast modulo a prime — without overflow, within time limits — is quite another. This post is the computational toolkit you'll reach for in almost every combinatorics contest problem.

The Binomial Theorem

Before diving into computation, let's state the theorem that gives binomial coefficients their name:

(x + y)ⁿ = Σₖ₌₀ⁿ C(n, k) · xᵏ · yⁿ⁻ᵏ

Setting x = y = 1 gives 2ⁿ = Σ C(n, k), confirming the row sum. Setting x = 1, y = −1 gives 0 = Σ (−1)ᵏ C(n, k), the alternating sum. These substitutions are a standard technique for evaluating sums involving binomial coefficients.

Useful corollaries

• Σ k·C(n,k) = n·2ⁿ⁻¹. Differentiate (1+x)ⁿ and set x = 1.
• Σ C(n,k)² = C(2n,n). Vandermonde's identity with m = n = r.
• Σ C(n,k)·2ᵏ = 3ⁿ. Set x = 2, y = 1.

Precomputing Factorials

The standard approach for computing C(n, r) mod p (where p is prime) requires two precomputed arrays of size MAXN:

• fac[i] = i! mod p
• ifac[i] = (i!)⁻¹ mod p (modular inverse of i!)

Then: C(n, r) = fac[n] · ifac[r] · ifac[n−r] mod p. Each query is O(1) after O(MAXN) precomputation.

const int MAXN = 2000005;
const long long MOD = 1e9 + 7;
long long fac[MAXN], ifac[MAXN];

long long pw(long long b, long long e, long long m = MOD) {
    long long r = 1; b %= m;
    for (; e > 0; e >>= 1) {
        if (e & 1) r = r * b % m;
        b = b * b % m;
    }
    return r;
}

void precompute() {
    fac[0] = 1;
    for (int i = 1; i < MAXN; i++)
        fac[i] = fac[i - 1] * i % MOD;

    // Single modular exponentiation, then propagate downward
    ifac[MAXN - 1] = pw(fac[MAXN - 1], MOD - 2);
    for (int i = MAXN - 2; i >= 0; i--)
        ifac[i] = ifac[i + 1] * (i + 1) % MOD;
}

long long C(int n, int r) {
    if (r < 0 || r > n) return 0;
    return fac[n] % MOD * ifac[r] % MOD * ifac[n - r] % MOD;
}

Complexity analysis

Modular Inverse for nCr mod p

The precomputation relies on computing modular inverses. For prime p, Fermat's little theorem gives us:

a⁻¹ ≡ aᵖ⁻² (mod p)

This works because aᵖ⁻¹ ≡ 1 (mod p) when gcd(a, p) = 1, so aᵖ⁻² · a ≡ 1 (mod p).

Alternative: linear inverse table

You can precompute inverses of 1 through N in O(N):

long long inv_table[MAXN];

void buildInverses() {
    inv_table[1] = 1;
    for (int i = 2; i < MAXN; i++)
        inv_table[i] = (MOD - MOD / i) * inv_table[MOD % i] % MOD;
}

This uses the identity: i⁻¹ ≡ −⌊p/i⌋ · (p mod i)⁻¹ (mod p). However, for computing binomial coefficients the factorial+inverse-factorial approach is more convenient.

Extended Euclidean (for composite moduli)

long long extgcd(long long a, long long b, long long &x, long long &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    long long x1, y1;
    long long g = extgcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

long long modinv(long long a, long long m) {
    long long x, y;
    long long g = extgcd(a, m, x, y);
    if (g != 1) return -1;  // inverse doesn't exist
    return (x % m + m) % m;
}

Lucas' Theorem

When n is very large (up to 10¹⁸) but p is small (a prime, say ≤ 10⁶), we can't precompute fac[n]. Instead, use Lucas' theorem:

C(n, r) mod p = Πᵢ C(nᵢ, rᵢ) mod p

where nᵢ and rᵢ are the digits of n and r in base p. If any rᵢ > nᵢ, the product is 0.

// Lucas' theorem implementation
// Precompute fac[] and ifac[] only up to p-1 (since ni, ri < p)
long long lucasC(long long n, long long r, long long p,
                 long long* fac_small, long long* ifac_small) {
    long long result = 1;
    while (n > 0 || r > 0) {
        int ni = n % p, ri = r % p;
        if (ri > ni) return 0;
        result = result % p * fac_small[ni] % p * ifac_small[ri] % p
                 * ifac_small[ni - ri] % p;
        n /= p;
        r /= p;
    }
    return result;
}

Limitations

• Only works when p is prime. For prime powers, use Andrew Granville's generalisation.
• If p is large (like 10⁹+7), there's no benefit — just use normal precomputation.

Large nCr with the Chinese Remainder Theorem

What if the modulus m is composite? Factor m = p₁^a₁ · p₂^a₂ · … · pₖ^aₖ, compute C(n, r) mod pᵢ^aᵢ for each prime power, then combine using CRT.

Computing C(n, r) mod pᵃ for a prime power is more involved — you need to extract the powers of p from the factorials. Here's the key idea:

Legendre's formula

The exponent of prime p in n! is:

v_p(n!) = ⌊n/p⌋ + ⌊n/p²⌋ + ⌊n/p³⌋ + …

// Exponent of prime p in n!
long long legendreExp(long long n, long long p) {
    long long e = 0;
    for (long long pk = p; pk <= n; pk *= p)
        e += n / pk;
    return e;
}

The exponent of p in C(n, r) is v_p(n!) − v_p(r!) − v_p((n−r)!). If this is ≥ 1, then C(n, r) ≡ 0 (mod p). Combined with Granville's method for the non-p part, you can compute C(n, r) mod pᵃ.

Implementation Template

Here is a battle-tested, copy-paste-ready template that covers the most common case: C(n, r) mod p where p is a large prime and n ≤ 2 × 10⁶.

#include <bits/stdc++.h>
using namespace std;

// ===================== Combinatorics Template =====================
namespace Combo {
    const long long MOD = 1e9 + 7;
    const int MAXN = 2000005;
    long long fac[MAXN], ifac[MAXN];

    long long pw(long long b, long long e, long long m = MOD) {
        long long r = 1; b %= m;
        for (; e > 0; e >>= 1) {
            if (e & 1) r = r * b % m;
            b = b * b % m;
        }
        return r;
    }

    void init() {
        fac[0] = 1;
        for (int i = 1; i < MAXN; i++)
            fac[i] = fac[i-1] * i % MOD;
        ifac[MAXN-1] = pw(fac[MAXN-1], MOD - 2);
        for (int i = MAXN-2; i >= 0; i--)
            ifac[i] = ifac[i+1] * (i+1) % MOD;
    }

    long long C(int n, int r) {
        if (r < 0 || r > n) return 0;
        return fac[n] % MOD * ifac[r] % MOD * ifac[n-r] % MOD;
    }

    long long P(int n, int r) {
        if (r < 0 || r > n) return 0;
        return fac[n] % MOD * ifac[n-r] % MOD;
    }

    long long catalan(int n) {
        return C(2*n, n) % MOD * pw(n+1, MOD-2) % MOD;
    }
}
// ==================================================================

// ===================== Lucas (small prime) ========================
namespace Lucas {
    long long fac_s[1000006], ifac_s[1000006];
    long long p_mod;

    void init(long long p) {
        p_mod = p;
        fac_s[0] = 1;
        for (int i = 1; i < p; i++)
            fac_s[i] = fac_s[i-1] * i % p;
        ifac_s[p-1] = Combo::pw(fac_s[p-1], p - 2, p);
        for (int i = (int)p - 2; i >= 0; i--)
            ifac_s[i] = ifac_s[i+1] * (i+1) % p;
    }

    long long C(long long n, long long r) {
        long long result = 1;
        while (n > 0 || r > 0) {
            int ni = n % p_mod, ri = r % p_mod;
            if (ri > ni) return 0;
            result = result * fac_s[ni] % p_mod
                     * ifac_s[ri] % p_mod
                     * ifac_s[ni - ri] % p_mod;
            n /= p_mod;
            r /= p_mod;
        }
        return result;
    }
}
// ==================================================================

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    Combo::init();

    int q; cin >> q;
    while (q--) {
        long long n, r;
        cin >> n >> r;
        // Use Combo::C for n <= 2e6 with MOD = 1e9+7
        // Use Lucas::C for very large n with small prime p
        cout << Combo::C(n, r) << "\n";
    }

    return 0;
}

Quick reference: which method to use

Practice Problems

1. CSES — Binomial Coefficients. Straightforward: answer q queries of C(n, r) mod 10⁹+7. Perfect for testing your precomputation template.
2. CSES — Bracket Sequences I. Count valid bracket sequences of length 2n → Catalan number = C(2n,n)/(n+1).
3. Codeforces 1359E — Modular Stability. Reduce to choosing k elements with a divisibility condition — uses C(n/d, k).
4. Codeforces 1696D — Permutation Graph. Though graph-themed, the counting part requires fast binomial coefficient queries.
5. AtCoder ABC 235 F — Variety of Digits. Digit DP combined with binomial coefficients and inclusion-exclusion.
6. Codeforces 1ившC (Div. 1) — Anu Has a Function. Bit-level analysis + counting with binomial coefficients.
7. Kattis — Binomial Coefficient (Large). Requires Lucas' theorem due to large n with small prime modulus.

Summary cheat sheet