Binomial coefficients

The binomial coefficient C(n, k) is the number of ways to choose $k$ items out of $n$ when the order doesn't matter. You compute it with Pascal's triangle when $n$ is at most a few thousand, and with factorials modulo a prime when $n$ is up to 10^6. To list the choices themselves rather than count them, use backtracking.

Explanation

Four friends, Ann, Ben, Cat and Dan, pick a team of 2 for a contest. Fig. 1 shows all the teams they can make: there are 6.

Fig. 1. There are 6 ways to choose 2 of 4 friends.

The first member can be any of the 4 friends and the second any of the 3 others, so there are 4 · 3 = 12 ordered picks. Each team appears among them exactly twice, as Ann then Ben and as Ben then Ann, so there are 12 / 2 = 6 teams. In general, choosing $k$ of $n$ in order gives n \cdot (n - 1) \cdot \ldots \cdot (n - k + 1) ordered picks, and each team appears once for each of its $k!$ orders. Dividing gives the formula

read "$n$ choose $k$". These numbers are called binomial coefficients.

A class of 30 picks a team of 15. The answer, C(30, 15), is 155117520, a nine-digit number. But the formula needs $30!$, which has 33 digits, and the largest standard integer type, long long, holds about 19.

The natural thing to try is to make the numbers smaller first: cancel $15!$, multiply only 16 \cdot 17 \cdot \ldots \cdot 30, then divide by $15!$. That product is about 2.0 \cdot 10^{20}, still too big: in long long it overflows, and dividing what is left by $15!$ gives −54279 instead of 155117520. Dividing as you go works better. Multiply by 16 and divide by 1, multiply by 17 and divide by 2, and so on; every step gives a whole number, and the result is right, but it still divides. Problems with large answers ask for them modulo a number such as 10^9 + 7, and the remainder of a quotient is not the quotient of the remainders. A method that only adds avoids the question.

To count with additions only, let Eve join the friends, so 5 friends choose a team of 2, and split the teams into groups instead of dividing. Every team either includes Ann or doesn't. If Ann is in, the other place goes to one of the 4 others: 4 teams. If Ann is out, both places go to the 4 others: 6 teams, as many as in Fig. 1. No team is in both groups, and none is missed, so there are 4 + 6 = 10 teams.

Fig. 2. Every team of 2 out of 5 friends either includes Ann or doesn't: 4 + 6 = 10.

The same split works for any number of friends and any team size. Pick one friend: the teams with that friend need k - 1 more members out of the other n - 1, and the teams without need all $k$ out of the other n - 1:

The splitting stops at the base cases: there is exactly 1 way to choose nobody and 1 way to choose everybody, and 0 ways to choose more friends than there are.

Write the counts in a table, one row for each number of friends $n$ and one column for each team size $k$. This table is Pascal's triangle, and filling it row by row is dynamic programming: each cell is computed once from two smaller ones.

The table is empty. We want row 5, choose 2.

Row 0: there is 1 way to choose nobody out of nobody.

Row 1: base cases, 1 way to choose nobody and 1 to choose everybody.

Row 2: base cases 1 and 1.

Row 2, choose 1: 1 + 1 = 2.

Row 3: base cases 1 and 1.

Row 3, choose 1: 1 + 2 = 3.

Row 3, choose 2: 2 + 1 = 3.

Row 4: base cases 1 and 1.

Row 4, choose 1: 1 + 3 = 4.

Row 4, choose 2: 3 + 3 = 6.

Row 4, choose 3: 3 + 1 = 4.

Row 5: base cases 1 and 1.

Row 5, choose 1: 1 + 4 = 5.

Row 5, choose 2: 4 + 6 = 10. There are 10 teams of 2 out of 5 friends.

Row 5, choose 3: 6 + 4 = 10.

Row 5, choose 4: 4 + 1 = 5. Row 5 is complete.

1/17

Fig. 3. Step through Pascal's triangle with the buttons or the arrow keys. Each new cell (amber) is the sum of the two blue cells above it.

The run computed every C(n, k) for $n$ up to 5, not only C(5, 2). The same loop up to row 30 gives C(30, 15) = 155117520. Every number added counts some of the teams, so that is the largest number it stores. Each row reads the same from both ends because choosing who goes is the same as choosing who stays: C(n, k) = C(n, n - k).

Large n modulo a prime

Many problems have $n$ up to 10^6 and ask for the answer modulo 10^9 + 7. The triangle would need about 5 \cdot 10^{11} cells, so we go back to the formula. If we take the remainder after every multiplication, the factorials stay small.

This brings back the division problem from the start. Take C(6, 2) = 15 modulo 13, which leaves 2. 6! = 720 leaves 5, and 2! \cdot 4! = 48 leaves 9. Dividing the remainders gives 5 / 9 = 0 in integer division, which is wrong.

Instead of dividing by 9, multiply by the modular inverse of 9: the number that gives remainder 1 when multiplied by 9. Modulo 13 it is 3, because 9 · 3 = 27 leaves 1. Then 5 · 3 = 15 leaves 2, the right answer. We store the inverse of every factorial: the inverse of $2!$ is 7 and of $4!$ is 6, so C(6, 2) leaves the same remainder as 5 · 7 · 6 = 210, which is 2. One inverse is enough to find them all. The inverse of $6!$ is 8, and each step down multiplies by the number being dropped: 8 · 6 leaves 9, 9 · 5 leaves 6, 6 · 4 leaves 11, 11 · 3 leaves 7. With the factorials and their inverses computed once, each C(n, k) is two multiplications.

Fig. 4. Modulo 13, C(6, 2) = 15 leaves 2: multiply the factorial of 6 by the inverse factorials of 2 and 4.

Pascal's triangle only adds, so any modulus works, but its n^2 cells limit $n$ to a few thousand: 5000 rows already take 200 MB. The factorial method reaches n = 10^6 and beyond, but it needs inverses, which exist only modulo a prime. Modulo 4, for example, 4! = 24 leaves 0, and no number multiplied by 0 leaves 1. So $4!$ has no inverse, and the factorial method gives 0 for C(4, 2) = 6, whose remainder is 2. The prime must also be larger than every $n$ you use. Modulo 13, $13!$ leaves 0, so with factorials up to 15 the inverse of $15!$ comes out as 0, every inverse factorial becomes 0, and even C(6, 2) comes out as 0. With 10^9 + 7 and n \le 10^6 this can't happen.

Implementation

The friends are the items, and a team is a choice of k of them. The table C is a 2D vector of long long, where C[n][k] sits in row n, column k.

1. State: C[n][k] is the number of ways to choose k of n items.

2. Base case: C[0][0] = 1 and C[n][0] = 1 for every n. Cells with k > n stay 0.

3. Order: rows n from 1 to N, inside each row k from 1 to n.

4. Transition: C[n][k] = (C[n - 1][k - 1] + C[n - 1][k]) % MOD.

5. Answer: C[n][k].

Since C starts filled with zeros, the last cell of a row needs no special case: C[n][n] adds C[n - 1][n - 1], which is 1, and C[n - 1][n], which is 0. Without % MOD and with N = 66, it gives exact values: row 66 is the last row where every value fits in long long. With N = 1000, C takes about 8 MB.

The input is n and k, with 0 \le k \le n \le 1000. For the five friends it is:

5 2

#include <iostream>
#include <vector>
using namespace std;

const long long MOD = 1000000007;
const int N = 1000;   // rows 0..N

int main() {
    // state: C[n][k] = the number of ways to choose k of n items; cells with k > n stay 0
    vector<vector<long long>> C(N + 1, vector<long long>(N + 1, 0));

    C[0][0] = 1;                                              // base case
    // order: rows n from 1 to N, inside a row k from 1 to n
    for (int n = 1; n <= N; n++) {
        C[n][0] = 1;                                          // base case: choose nobody
        for (int k = 1; k <= n; k++)
            C[n][k] = (C[n - 1][k - 1] + C[n - 1][k]) % MOD;  // transition; remove % MOD and set N = 66 for exact values
    }

    int n, k;
    cin >> n >> k;
    cout << C[n][k] << "\n";                                  // answer
}

There are O(N^2) cells, and each takes one addition, so the time and the memory grow as the square of N, O(N^2).

Factorials modulo a prime

The program computes fact and inv_fact once, up to MAXN.

1. Factorials: fact[0] = 1 and fact[i] = fact[i - 1] * i % MOD.

2. Inverse factorials: inv_fact[MAXN] = power(fact[MAXN], MOD - 2, MOD), then downwards inv_fact[i - 1] = inv_fact[i] * i % MOD. For a prime MOD larger than MAXN, x to the power MOD - 2 is the inverse of x (Fermat's little theorem), and the inverse of (i - 1)! is i times the inverse of i!. The modular inverse article explains why.

3. Combine: C(n, k) returns 0 if k < 0 or k > n, and fact[n] * inv_fact[k] % MOD * inv_fact[n - k] % MOD otherwise.

The function power is the one from binary exponentiation. The input is the number of queries q, then q pairs n k with n \le 10^6:

3
6 2
30 15
1000000 500000

#include <iostream>
#include <vector>
using namespace std;

const long long MOD = 1000000007;   // a prime
const int MAXN = 1000000;           // the largest n in C(n, k)

long long power(long long x, long long n, long long m) {
    long long result = 1 % m;
    x %= m;
    while (n > 0) {
        if (n % 2 == 1) result = result * x % m;
        x = x * x % m;
        n /= 2;
    }
    return result;
}

vector<long long> fact(MAXN + 1), inv_fact(MAXN + 1);

long long C(int n, int k) {
    if (k < 0 || k > n) return 0;                                 // no way to choose
    return fact[n] * inv_fact[k] % MOD * inv_fact[n - k] % MOD;   // combine
}

int main() {
    // factorials
    fact[0] = 1;
    for (int i = 1; i <= MAXN; i++)
        fact[i] = fact[i - 1] * i % MOD;

    // inverse factorials: one power, then downwards
    inv_fact[MAXN] = power(fact[MAXN], MOD - 2, MOD);
    for (int i = MAXN; i > 0; i--)
        inv_fact[i - 1] = inv_fact[i] * i % MOD;

    int q;
    cin >> q;
    while (q--) {
        int n, k;
        cin >> n >> k;
        cout << C(n, k) << "\n";
    }
}

The precomputation takes O(\text{MAXN} + \log \text{MOD}) time and $O(\text{MAXN})$ memory, and each query takes $O(1)$.

Common mistakes

• Filling the triangle in int. C(34, 17) = 2333606220 doesn't fit, so use long long.

• Multiplying three remainders before reducing, as in fact[n] * inv_fact[k] * inv_fact[n - k] % MOD. Each factor is below 10^9 + 7, so the product of two fits in long long but the product of three doesn't: with this line the program prints a negative number even for C(6, 2). Take % MOD after each multiplication.

• Calling C(n, k) without checking that 0 \le k \le n. When n - k comes from a formula, as in counting paths, it can be negative: C(3, 5) reads inv_fact[-2], so return 0 first.

• Sizing fact and inv_fact by the bound on the input instead of the largest first argument. For C(n + m, n) with n, m \le 10^6 they need 2 \cdot 10^6 + 1 entries, or fact[n + m] reads outside them.

• Keeping one row of the triangle and updating row[k] += row[k - 1] with k going up. Row 1 2 2 comes out for n = 2, so C(2, 2) = 2, because row[1] was already overwritten. Go from k = n down to 1.

Similar algorithms

• Grid paths fill a table by additions too. On an empty grid the number of right and down paths is a binomial coefficient, C(n + m - 2, n - 1); with blocked cells you need the table.

• Stars and bars counts choices with repetition, such as identical items put into boxes. The answer is a binomial coefficient of a changed problem, C(n + k - 1, k - 1).

• Catalan numbers count correct bracket sequences and similar structures. They equal C(2n, n) / (n + 1) and are computed with the same factorials.

• Backtracking lists the choices one by one.

• Big integers give exact values beyond row 66.

Practice