Factorials and binomials modulo p

Why this matters when you build with AI

Many contest problems ask for a count, and the count is a binomial coefficient C(n, k). It is the number of ways to choose k items from n items when the order does not matter. The statement usually says "print the answer modulo 998244353", because the true number is far too large to store. An AI tool asked for this often writes a loop that multiplies the numbers from 1 to n and then divides. That code is correct for small n and wrong for some large n, and the wrong cases are easy to miss.

Here is a failure from our runs. We wrote the factorial version with no guard, and asked it for C(8, 1) modulo 7. The true value is 8, and 8 leaves remainder 1 when divided by 7. The program printed 0. Nothing crashed, and 0 looks like a normal answer. Sample tests often use n below the modulus, so the bug stays hidden until a large test.

This lesson covers the standard method, the two limits that keep it correct, one worked problem, and a test you can run on your own machine. The code is C++17. We compiled it with g++ -std=c++17 -O2 -Wall -Wextra using GCC 15.2. The output lines below are from our runs. Your output should match them exactly, because the arithmetic does not depend on the machine.

Definitions

C(n, k) is the binomial coefficient. For k from 0 to n, it counts the k-element subsets of an n-element set. For example, C(5, 2) is 10, because the five items A, B, C, D, E have 10 pairs.

A modulus is the number we divide by when we want the remainder. "x modulo p" means the remainder of x divided by p. A prime is a number whose only divisors are 1 and itself. The primes used most in contests are 998244353 and 1000000007. The number 6 is not prime, because 2 and 3 divide it.

The factorial n! is 1 times 2 times ... times n, and 0! is 1. So 5! is 120 and 2! is 2.

The formula and the idea

The binomial coefficient is C(n, k) = n! / (k! times (n - k)!). With n = 5 and k = 2, that is 120 / (2 times 6), which is 10.

Modulo a prime p, division becomes multiplication by an inverse. The inverse of x is a number y such that x times y leaves remainder 1 modulo p. For a prime p, and any x that p does not divide, Fermat's little theorem gives the inverse as x raised to the power p - 2, taken modulo p. The module on number theory has a full lesson on this, named "Modular inverse and division".

So we need two tables. The first is fact[x] = x! modulo p. The second is inv_fact[x], the inverse of x! modulo p. Then:

C(n, k) modulo p = fact[n] times inv_fact[k] times inv_fact[n - k], taken modulo p.

Each table entry is found from the one before it. We compute inv_fact[N] once with one exponentiation. Then we go down: inv_fact[x - 1] = inv_fact[x] times x, modulo p. The reason is that x! = x times (x - 1)!, so the inverse of (x - 1)! is the inverse of x! times x.

Now a worked example with p = 13 and C(6, 2). Work each factorial modulo 13.

fact[6] = 720. Since 13 times 55 is 715, fact[6] leaves 5.
fact[2] = 2 and fact[4] = 24, which leaves 11.
inv_fact[2] = 7, because 2 times 7 is 14, which leaves 1.
inv_fact[4] = 6, because 11 times 6 is 66, which leaves 1.

Then C(6, 2) modulo 13 is 5 times 7 times 6, which is 210. Since 13 times 16 is 208, 210 leaves 2. Check it directly: C(6, 2) is 15, and 15 leaves 2 modulo 13. The two answers agree.

The inverse of x! exists only when x! is not divisible by p. A product of numbers is divisible by a prime p exactly when one of its factors is. So x! is safe when every number from 1 to x is below p. This gives the first rule of this lesson: the table size N must be less than p.

The code

The struct Binom builds both tables in O(N) time and answers each query in O(1) time. The main function prints a small Pascal table and two traps. Copy it, compile it, and run it.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
 
ll pow_mod(ll a, ll e, ll m) {
    ll r = 1 % m;
    a %= m;
    while (e > 0) {
        if (e & 1) r = r * a % m;
        a = a * a % m;
        e >>= 1;
    }
    return r;
}
 
// Factorials and inverse factorials modulo a prime p, for 0 <= n <= N.
// Needs N < p, so that no factorial up to N is divisible by p.
struct Binom {
    ll p;
    int N;
    vector<ll> fact, inv_fact;
    Binom(int N_, ll p_) : p(p_), N(N_), fact(N_ + 1), inv_fact(N_ + 1) {
        if (N >= p) throw invalid_argument("need N < p");
        fact[0] = 1;
        for (int i = 1; i <= N; i++) fact[i] = fact[i - 1] * i % p;
        inv_fact[N] = pow_mod(fact[N], p - 2, p);
        for (int i = N; i >= 1; i--) inv_fact[i - 1] = inv_fact[i] * i % p;
    }
    // C(n, k) mod p. Returns 0 when k < 0 or k > n.
    ll C(ll n, ll k) const {
        if (k < 0 || k > n) return 0;
        if (n > N) throw out_of_range("n above the table");
        return fact[n] * inv_fact[k] % p * inv_fact[n - k] % p;
    }
};
 
// The trap: build the same tables without the N < p check.
// Shows what goes wrong when N >= p.
ll naive_C(int N, ll p, ll n, ll k) {
    vector<ll> fact(N + 1), inv(N + 1);
    fact[0] = 1;
    for (int i = 1; i <= N; i++) fact[i] = fact[i - 1] * i % p;
    inv[N] = pow_mod(fact[N], p - 2, p);   // fact[N] is 0 when p <= N, so this is 0
    for (int i = N; i >= 1; i--) inv[i - 1] = inv[i] * i % p;
    if (k < 0 || k > n) return 0;
    return fact[n] * inv[k] % p * inv[n - k] % p;
}
 
int main() {
    const ll p = 13;
    Binom B(12, p);
    cout << "Pascal rows 0..6 from the factorial formula, mod 13:\n";
    for (int n = 0; n <= 6; n++) {
        for (int k = 0; k <= n; k++) cout << B.C(n, k) << (k == n ? '\n' : ' ');
    }
    cout << "C(12, 6) mod 13 = " << B.C(12, 6) << "\n";
    cout << "C(5, 7) mod 13 = " << B.C(5, 7) << " (k > n)\n";
    cout << "C(5, -1) mod 13 = " << B.C(5, -1) << " (k < 0)\n";
 
    cout << "\nThe trap, p = 7, table up to 8 (N >= p):\n";
    cout << "naive C(8, 1) mod 7 = " << naive_C(8, 7, 8, 1) << "   (true value 8 mod 7 = 1)\n";
    cout << "naive C(7, 0) mod 7 = " << naive_C(8, 7, 7, 0) << "   (true value 1)\n";
    try {
        Binom bad(8, 7);
    } catch (const invalid_argument& e) {
        cout << "Binom(8, 7) throws: " << e.what() << "\n";
    }
    return 0;
}

Compiled and run with g++ -std=c++17 -O2 -Wall -Wextra, it printed:

Pascal rows 0..6 from the factorial formula, mod 13:
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
1 6 2 7 2 6 1
C(12, 6) mod 13 = 1
C(5, 7) mod 13 = 0 (k > n)
C(5, -1) mod 13 = 0 (k < 0)
 
The trap, p = 7, table up to 8 (N >= p):
naive C(8, 1) mod 7 = 0   (true value 8 mod 7 = 1)
naive C(7, 0) mod 7 = 0   (true value 1)
Binom(8, 7) throws: need N < p

Check the row 6 values by hand. C(6, 2) is 15, which leaves 2 modulo 13. C(6, 3) is 20, which leaves 7. Both match the printed row. C(12, 6) is 924, and 13 times 71 is 923, so the remainder is 1. That also matches.

The n >= p trap

The naive version failed because it built the table up to N = 8 with p = 7. The number 7 divides 7!, so fact[7] and fact[8] are both 0 modulo 7. The inverse of 0 is found as pow(0, 5, 7), which is 0. Every table value from index 7 upward is 0. Any product that uses one of them is 0. So C(8, 1) printed 0, and the true remainder is 1.

The same thing happens for C(7, 0). The true value is 1, but the table entry fact[7] is 0, so the answer is 0.

A fix that returns 0 whenever n >= p is also wrong. For example, C(8, 1) is 8, which is not 0 modulo 7. For n >= p, the next lesson gives the correct method, Lucas' theorem, which is named "Pascal's triangle and Lucas' theorem" in this module. The Binom struct refuses to build when N >= p, and its error message is the first thing to look for when a table answer is strange.

Out-of-range k

A binomial coefficient is 0 when k is below 0 or above n, because there are no such subsets. The check if (k < 0 || k > n) return 0; must come before any array index. A negative k would read inv_fact[-1], which is outside the array. C++ does not check array bounds, so the program may print garbage or crash.

The struct also guards n. If n is above the table size N, it throws out_of_range. This is better than reading past the end of the vector.

The same ideas in Python

Python has pow(x, -1, p), which returns the inverse of x modulo p. It needs Python 3.8 or later. The version below builds the same tables for p = 13 and N = 12. It does the same linear work as the C++ version, but Python is slower per step, so use the C++ version for a contest with large N.

Python · editable

Its output is:

[1, 6, 2, 7, 2, 6, 1]
1 1

Python's comb is exact, so the second line checks the first against an independent value.

Worked problem: answer many queries

Statement. You get a number of queries Q. Each query gives two integers n and k. Print C(n, k) modulo 998244353. Print 0 when k is negative or larger than n.

Limits: 1 to 100000 queries, 0 <= n <= 200000, and -10 <= k <= 200010.

Sample input:

5
5 2
10 0
4 7
200000 100000
3 -1

Sample output:

Solution. The modulus 998244353 is prime, and every n is at most 200000, which is below the modulus. So the table in Binom is valid. Build it once with N = 200000. Then answer each query with the check for k and one formula. The total time is O(N + Q), and the memory is two arrays of 200001 64-bit numbers, about 3.2 MB.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
 
ll pow_mod(ll a, ll e, ll m) {
    ll r = 1 % m;
    a %= m;
    while (e > 0) {
        if (e & 1) r = r * a % m;
        a = a * a % m;
        e >>= 1;
    }
    return r;
}
 
struct Binom {
    ll p;
    int N;
    vector<ll> fact, inv_fact;
    Binom(int N_, ll p_) : p(p_), N(N_), fact(N_ + 1), inv_fact(N_ + 1) {
        if (N >= p) throw invalid_argument("need N < p");
        fact[0] = 1;
        for (int i = 1; i <= N; i++) fact[i] = fact[i - 1] * i % p;
        inv_fact[N] = pow_mod(fact[N], p - 2, p);
        for (int i = N; i >= 1; i--) inv_fact[i - 1] = inv_fact[i] * i % p;
    }
    ll C(ll n, ll k) const {
        if (k < 0 || k > n) return 0;
        if (n > N) throw out_of_range("n above the table");
        return fact[n] * inv_fact[k] % p * inv_fact[n - k] % p;
    }
};
 
int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    const ll MOD = 998244353;
    const int MAXN = 200000;
    Binom B(MAXN, MOD);
    int q;
    cin >> q;
    while (q--) {
        ll n, k;
        cin >> n >> k;
        cout << B.C(n, k) << "\n";
    }
    return 0;
}

Compiled and run on the sample input, it printed:

Each line matches the sample. The third query has k = 7 above n = 4, so it is 0. The fifth query has k = -1, so it is 0. To check the large value, we ran math.comb(200000, 100000) % 998244353 in Python, which is exact big-integer arithmetic. It printed 935629711, the same as the program.

Testing against brute force

A brute force counts the subsets directly. It walks every bitmask of n bits and counts the masks with k bits set. This is slow, so it is only for n up to about 20.

// Count subsets of size k of {0..n-1} by trying every bitmask. Only for n <= 20.
ll brute_count(int n, int k) {
    ll cnt = 0;
    for (int mask = 0; mask < (1 << n); mask++)
        if (__builtin_popcount(mask) == k) cnt++;
    return cnt;
}

Because the brute force is exact for small n, it checks the table at every (n, k) pair where it can. We ran the test on the moduli 13, 17, 101 and 1000000007, for every n up to the smaller of p - 1 and 18, and every k from -2 to n + 2. It agreed on all 896 cases. We also compared the table against an independent Pascal triangle with the rule C(n, k) = C(n - 1, k - 1) + C(n - 1, k), taken modulo 1000000007 for n up to 500. It agreed on 2000 random pairs.

Complexity

Precomputing the tables takes O(N) multiplications plus one exponentiation, which is O(N + log p). Each query takes O(1) time. The memory is O(N). Here N is the largest n that the problem allows, and it must be below p.

Common mistakes

Using a table with N >= p. The symptom is a 0 where the true answer is not 0, such as the naive C(8, 1) mod 7 printing 0 instead of 1. The fix is to check N < p, or to use Lucas' theorem when n is large.
Computing the inverse with the same exponent for a non-prime modulus. Fermat's rule needs a prime p. For p = 6, the code pow(2, 6 - 2, 6) returns 4, and 2 times 4 is 8, which leaves 2, not 1. So 4 is not the inverse of 2. Python's built-in gives a clear error, as in pow(2, -1, 6), which raises ValueError with the message "base is not invertible for the given modulus". C++ does not raise an error, so the wrong number goes on silently.
Reading inv_fact[k] before checking k. A negative k or a k above the table reads outside the array. C++ does not check this, so the program may print a wrong value or crash.
Multiplying two values without reducing modulo p first, or storing them in int. Two values near 998244353 multiply to about 10^18. That fits in a 64-bit long long, which holds up to about 9.2 times 10^18. It does not fit in a 32-bit int, and an overflow in C++ is undefined behaviour, so the result is not reliable.
Rebuilding the tables for every query. The tables take O(N) time, so rebuilding them Q times costs O(N times Q). Build them once.

A prompt to check AI output

When an AI tool writes a binomial table for you, paste its code into a new chat with this prompt.

Here is my C++ code for C(n, k) modulo a prime p. Check it against these four rules and mark each one pass or fail.
1. The table size N is at most p - 1.
2. k < 0 and k > n return 0 before any array index is used.
3. Every product is reduced modulo p, and the variables are 64-bit.
4. The inverse of N! comes from pow(fact[N], p - 2, p), and the other inverses come from going down the table.
For each fail, give the smallest input that breaks it, and the value a brute-force count gives for that input. Do not change the code yet.

A good answer marks each rule and gives a counterexample for each fail. For a table with N = 8 and p = 7, the counterexample is C(8, 1), with true value 1. Check the answer by running the counterexample yourself, because the model can state a wrong value with confidence.

Checklist

Is the modulus prime? If not, the factorial method does not apply.
Is N below p? If n can reach p, use Lucas' theorem instead.
Does the code return 0 for k < 0 and for k > n before any array index?
Is every product reduced modulo p, and is it in a 64-bit type?
Is the table built once, outside the query loop?
Did you compare small n with a brute-force count?

Practice

Each practice problem has a hint. Try it before you read the answer under the line "Show the answer".

Compute C(10, 4) modulo 11 using the factorial method, by hand. Hint: first find C(10, 4) as an integer, then reduce it.

Show the answer C(10, 4) is 210. Since 11 times 19 is 209, 210 leaves 1 modulo 11. Python's comb(10, 4) % 11 also gives 1.

Why does Binom(12, 13) build but Binom(13, 13) throw? Hint: which number between 1 and 13 is divisible by 13?

Show the answer The number 13 is a factor of 13!, so 13! is 0 modulo 13. The table up to N = 13 would need the inverse of 13!, which does not exist. For p = 13 the largest allowed table size is N = 12.

Compute C(100, 50) modulo 997 with Binom and N = 100. Hint: 100 is below 997, so the table is valid.

Show the answer The answer is 28. Python's comb(100, 50) % 997 also gives 28. To check the program, change the constants in the sample code to N = 100 and p = 997, and print B.C(100, 50).

A colleague builds the inverse of every x from 1 to N with its own call to pow_mod(x, p - 2, p). Is the answer wrong, and what is the cost? Hint: compare the number of exponentiations with the number of table entries.

Show the answer The answers are correct, because each pow_mod call gives the inverse of x when p is prime. The cost is about N exponentiations, and each one takes about 30 multiplications for p near 10^9. The table method needs one exponentiation and N multiplications in total.

Exercise: write your own table and check it

Write a function that builds the factorial and inverse factorial tables for a prime p and a table size N < p. Then write a brute-force count for small n, and check the table against it for every n from 0 to 16 and every k from -2 to n + 2. Use p = 13, 17, 101 and 1000000007. The check passes when every pair agrees. Keep the brute-force code and the table code separate, so that the brute force does not depend on the table.

Show the answer

The full test program is below. It uses the same Binom struct, so it checks the table you would submit. It also has a random check against a Pascal triangle. We compiled it with g++ -std=c++17 -O2 -Wall -Wextra and ran it. Its output is below the program.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
 
ll pow_mod(ll a, ll e, ll m) {
    ll r = 1 % m;
    a %= m;
    while (e > 0) {
        if (e & 1) r = r * a % m;
        a = a * a % m;
        e >>= 1;
    }
    return r;
}
 
struct Binom {
    ll p;
    int N;
    vector<ll> fact, inv_fact;
    Binom(int N_, ll p_) : p(p_), N(N_), fact(N_ + 1), inv_fact(N_ + 1) {
        if (N >= p) throw invalid_argument("need N < p");
        fact[0] = 1;
        for (int i = 1; i <= N; i++) fact[i] = fact[i - 1] * i % p;
        inv_fact[N] = pow_mod(fact[N], p - 2, p);
        for (int i = N; i >= 1; i--) inv_fact[i - 1] = inv_fact[i] * i % p;
    }
    ll C(ll n, ll k) const {
        if (k < 0 || k > n) return 0;
        if (n > N) throw out_of_range("n above the table");
        return fact[n] * inv_fact[k] % p * inv_fact[n - k] % p;
    }
};
 
// Count subsets of size k of {0..n-1} by trying every bitmask. Only for n <= 20.
ll brute_count(int n, int k) {
    ll cnt = 0;
    for (int mask = 0; mask < (1 << n); mask++)
        if (__builtin_popcount(mask) == k) cnt++;
    return cnt;
}
 
int main() {
    mt19937 rng(12345);
    long long checked = 0;
    for (ll p : {13LL, 17LL, 101LL, 1000000007LL}) {
        int N = (p < 1000000007LL) ? (int)min<ll>(p - 1, 18) : 18;
        Binom B(N, p);
        for (int n = 0; n <= N; n++) {
            for (int k = -2; k <= n + 2; k++) {
                ll want = (k < 0 || k > n) ? 0 : brute_count(n, k) % p;
                ll got = B.C(n, k);
                if (want != got) {
                    cout << "MISMATCH p=" << p << " n=" << n << " k=" << k << " want=" << want << " got=" << got << "\n";
                    return 1;
                }
                checked++;
            }
        }
    }
    cout << "brute force agrees on " << checked << " (p, n, k) cases\n";
 
    for (int trial = 0; trial < 2000; trial++) {
        ll p = 1000000007LL;
        int N = 500;
        Binom B(N, p);
        static vector<vector<ll>> P;
        if (P.empty()) {
            P.assign(N + 1, vector<ll>(N + 1, 0));
            for (int n = 0; n <= N; n++) {
                P[n][0] = 1;
                for (int k = 1; k <= n; k++) P[n][k] = (P[n - 1][k - 1] + (k <= n - 1 ? P[n - 1][k] : 0)) % p;
            }
        }
        int n = rng() % (N + 1);
        int k = (int)(rng() % (n + 3)) - 1;
        ll want = (k < 0 || k > n) ? 0 : P[n][k];
        if (want != B.C(n, k)) {
            cout << "RANDOM MISMATCH n=" << n << " k=" << k << "\n";
            return 1;
        }
    }
    cout << "2000 random cases agree with Pascal DP mod 1000000007\n";
    return 0;
}

brute force agrees on 896 (p, n, k) cases
2000 random cases agree with Pascal DP mod 1000000007

If your table disagrees, the first MISMATCH line gives the smallest failing n and k. Start with that pair, and look for the N < p rule and the k check.

Summary

For a prime p and a table size below p, build the factorial and inverse factorial tables once. Then answer each query in constant time, after the k check. Above that size, the table method fails, and the next lesson gives the method for large n.