Tổ hợp và đếm cơ bản

Rất nhiều bài lập trình thi đấu hỏi "có bao nhiêu cách...". Duyệt trâu từng cấu hình tốn thời gian tỉ lệ với chính số cách (thường là $2^n$, $n!$ hay lớn hơn), trong khi một công thức tổ hợp cho đáp án trong $O(1)$ sau khi tiền xử lý $O(n)$. Trang này đi từ hai quy tắc đếm cơ bản đến hoán vị, chỉnh hợp, tổ hợp, cách tính $\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ theo modulo nguyên tố và modulo hợp số, bài toán chia kẹo và số Catalan.

Ý tưởng / Trực giác

Hai quy tắc nền tảng

• Quy tắc cộng: nếu các cấu hình chia được thành những nhóm rời nhau (không cấu hình nào thuộc hai nhóm), số cấu hình bằng tổng số cấu hình của từng nhóm. Đây là cơ sở của mọi lời giải "xét từng trường hợp".
• Quy tắc nhân: nếu một cấu hình được dựng qua nhiều bước liên tiếp, và số lựa chọn ở mỗi bước không phụ thuộc vào lựa chọn cụ thể ở các bước trước, thì số cấu hình bằng tích số lựa chọn. Ví dụ: xâu độ dài $n$ trên bảng chữ cái $c$ ký tự có $c^n$ cách.

Hai điều kiện in đậm là chỗ hay sai nhất. Cộng hai nhóm chồng lên nhau thì đếm trùng (khi đó cần bao hàm loại trừ); nhân khi số lựa chọn thay đổi theo bước trước thì sai hoàn toàn.

Một mẹo đi kèm là đếm phần bù: đếm cấu hình "tốt" khó thì đếm tất cả trừ đi cấu hình "xấu".

Hoán vị, chỉnh hợp, tổ hợp

Chọn và xếp $k$ phần tử từ $n$ phần tử phân biệt:

Dòng thứ ba là mẫu lập luận quan trọng nhất của cả trang: đếm có thứ tự rồi chia cho số lần mỗi cấu hình bị đếm lặp. Nó chỉ đúng khi mọi cấu hình bị đếm lặp cùng một số lần.

Tam giác Pascal: tổ hợp chỉ bằng phép cộng

Xét phần tử cuối cùng: một tập $k$ phần tử của $\{1..n\}$ hoặc không chứa $n$ (chọn đủ $k$ từ $n-1$ phần tử còn lại) hoặc chứa $n$ (chọn thêm $k-1$). Theo quy tắc cộng:

$$\left(\genfrac{}{}{0ex}{}{n}{k}\right)=\left(\genfrac{}{}{0ex}{}{n-1}{k-1}\right)+\left(\genfrac{}{}{0ex}{}{n-1}{k}\right),\left(\genfrac{}{}{0ex}{}{n}{0}\right)=\left(\genfrac{}{}{0ex}{}{n}{n}\right)=1.$$

Công thức này không có phép chia, nên đúng với mọi modulo, kể cả modulo hợp số hay khi đếm chính xác bằng số lớn. Cái giá là $O(n^2)$ phép cộng, chỉ dùng được khi $n$ cỡ vài nghìn. Bộ nhớ cũng cần để ý: bảng đầy đủ long long C[2001][2001] đã chiếm khoảng $32$ MB. Nếu chỉ cần một hàng $n$ thì giữ một mảng một chiều và cập nhật từ phải sang trái, bộ nhớ còn $O(n)$.

Vì sao $\left(\genfrac{}{}{0ex}{}{n}{k}\right)modp$ cần nghịch đảo, và vì sao là Fermat

Khi $n$ tới ${10}^6$ mà cần nhiều truy vấn, ta dùng $\left(\genfrac{}{}{0ex}{}{n}{k}\right)=n!·(k!{)}^{-1}·((n-k)!{)}^{-1}$. Trong số học đồng dư không có phép chia, chỉ có nhân với nghịch đảo: $x^{-1}$ là số thỏa $x·x^{-1}\equiv 1\left(\mathrm{mod}p\right)$.

Với $p$ nguyên tố và $x\not\equiv 0$, định lý Fermat nhỏ cho $x^{p-1}\equiv 1$, nên

$$x^{-1}\equiv x^{p-2}\left(\mathrm{mod}p\right),$$

tính bằng lũy thừa nhanh $O(\log p)$. Nghịch đảo tồn tại khi và chỉ khi $gcd(x,p)=1$. Mọi thừa số trong $k!$ đều $<p$ khi $k<p$, nên $k!$ không chia hết cho $p$ và luôn nghịch đảo được. Các modulo quen thuộc ${10}^9+7$ và $998244353$ được chọn chính vì chúng nguyên tố và lớn hơn mọi $n$ trong đề.

Mẹo tiết kiệm: chỉ cần một lần lũy thừa. Tính inv_fact[N] $=(N!{)}^{p-2}$, rồi đi lùi: ${\displaystyle \frac{1}{(i-1)!}}={\displaystyle \frac{1}{i!}}·i$.

Chia kẹo (stars and bars)

Đếm số nghiệm nguyên không âm của $x_1+x_2+\cdots +x_k=n$ (chia $n$ viên kẹo giống nhau cho $k$ em). Xếp $n$ ngôi sao * và $k-1$ vách | thành một hàng: vách chia hàng sao thành $k$ nhóm, số sao nhóm thứ $i$ chính là $x_i$. Đây là song ánh, nên đáp án là số cách chọn vị trí cho $k-1$ vách trong $n+k-1$ ô:

$$\left(\genfrac{}{}{0ex}{}{n+k-1}{k-1}\right).$$

Nếu yêu cầu $x_i\ge 1$ thì phát trước mỗi em một viên, còn $n-k$ viên chia tự do: $\left(\genfrac{}{}{0ex}{}{n-1}{k-1}\right)$. Tổng quát, $x_i\ge a_i$ thì trừ trước $\sum a_i$. Ràng buộc chặn trên $x_i\le b_i$ thì cần bao hàm loại trừ.

Một biến thể hay gặp: số dãy không giảm độ dài $n$ với giá trị trong $[1,K]$ bằng $\left(\genfrac{}{}{0ex}{}{n+K-1}{n}\right)$ (đặt $x_v$ là số lần giá trị $v$ xuất hiện).

Số Catalan và nguyên lý phản xạ

Số dãy ngoặc đúng gồm $n$ cặp ngoặc là số Catalan

$$C_n=\frac{1}{n+1}\left(\genfrac{}{}{0ex}{}{2n}{n}\right)=\left(\genfrac{}{}{0ex}{}{2n}{n}\right)-\left(\genfrac{}{}{0ex}{}{2n}{n+1}\right):1,1,2,5,14,42,132,429,\dots$$

Coi ( là bước lên $+1$, ) là bước xuống $-1$. Dãy ngoặc đúng là đường đi $2n$ bước từ độ cao $0$ về $0$ không bao giờ xuống dưới $0$. Tổng số đường từ $0$ về $0$ (bất kể điều kiện) là $\left(\genfrac{}{}{0ex}{}{2n}{n}\right)$. Đường xấu là đường chạm độ cao $-1$. Lấy điểm đầu tiên nó chạm $-1$, lật đối xứng toàn bộ phần phía trước qua đường $y=-1$: đường mới xuất phát từ $-2$ và kết thúc ở $0$, tức có $n+1$ bước lên và $n-1$ bước xuống. Phép lật này là song ánh giữa đường xấu và mọi đường từ $-2$ tới $0$, nên số đường xấu là $\left(\genfrac{}{}{0ex}{}{2n}{n+1}\right)$. Trừ đi được công thức trên.

Nguyên lý phản xạ dùng được cả khi đã biết trước một tiền tố: nếu còn $r$ ký tự phải điền, độ cao hiện tại là $b\ge 0$, cần $d=\frac{r+b}{2}$ dấu ) thì đáp án là $\left(\genfrac{}{}{0ex}{}{r}{d}\right)-\left(\genfrac{}{}{0ex}{}{r}{d+1}\right)$.

Các bài quy về Catalan (đáp án đều là $C_n$):

• số dãy ngoặc đúng $n$ cặp; số đường đi trên lưới $n\times n$ từ góc dưới trái lên góc trên phải, chỉ đi phải/lên, không vượt lên trên đường chéo;
• số cây nhị phân có $n$ nút; số cách đặt ngoặc đầy đủ một tích $n+1$ thừa số (biểu thức $k$ toán hạng có $C_{k-1}$ cách);
• số cách chia đa giác lồi $n+2$ đỉnh thành tam giác bằng các đường chéo không cắt nhau;
• số hoán vị $n$ phần tử tránh một mẫu độ dài 3 cố định bất kỳ (ví dụ mẫu 321: không có 3 phần tử giảm dần; hay mẫu 231: sắp xếp được bằng một ngăn xếp);
• công thức truy hồi $C_0=1$, $C_{n+1}={\sum }_{i=0}^n{C}_i{C}_{n-i}$ (tách theo vị trí dấu ) khớp với dấu ( đầu tiên), không có phép chia nên dùng được với mọi modulo khi $n$ nhỏ.

Dấu hiệu nhận ra Catalan: tính tay đáp án cho $n=0,1,2,3,4$ ra $1,1,2,5,14$.

Khi modulo là hợp số

Với modulo như ${10}^9$, $2^{32}$ hay $m$ cho trong input, cách giai thừa + Fermat sai: $x^{m-2}$ không phải nghịch đảo, và tệ hơn, $k!$ có thể chia hết cho một ước nguyên tố của $m$ nên không có nghịch đảo nào. Ví dụ với $m=12$: $4!=24\equiv 0$. Có hai hướng xử lý:

1. Phân tích đáp án ra thừa số nguyên tố (công thức Legendre). Số mũ của nguyên tố $p$ trong $x!$ là $$v_p(x!)=\lfloor \frac{x}{p}\rfloor +\lfloor \frac{x}{p^2}\rfloor +\lfloor \frac{x}{p^3}\rfloor +\cdots$$ (có $\lfloor x/p\rfloor$ số chia hết cho $p$, trong đó $\lfloor x/p^2\rfloor$ số góp thêm một thừa số nữa, ...). Vậy $\left(\genfrac{}{}{0ex}{}{n}{k}\right)={\prod }_{p\le n}{p}^{e_p}$ với $e_p=v_p(n!)-v_p(k!)-v_p((n-k)!)$. Chỉ còn phép nhân, không còn phép chia. Hợp với một vài giá trị cần tính, $n$ tới cỡ ${10}^7$ (cần sàng nguyên tố).
2. Tách các ước nguyên tố của modulo ra riêng. Phân tích $m=p_1^{a_1}\cdots p_t^{a_t}$ (với $m\le {10}^9$ thì $t\le 9$). Viết mỗi số $i$ thành (phần gồm các $p_j$) $\times$ (phần nguyên tố cùng nhau với $m$). Phần sau nghịch đảo được bằng Euclid mở rộng; phần trước chỉ lưu số mũ rồi cộng trừ. Hợp với nhiều truy vấn.

Với modulo lớn cỡ $2^{32}$, tích hai số dư đã vượt $2^{63}$ nên phép a * b % m bằng long long bị tràn. Hai đoạn code ở mục Cài đặt nhân qua __int128 nên phép nhân đúng với mọi $m<2^{62}$, kể cả $m=2^{32}$. Riêng hướng 2 còn phải phân tích $m$: code mẫu chia thử tới $\sqrt{m}$, chỉ nhanh khi $m$ tới khoảng ${10}^{12}$ (cỡ ${10}^6$ bước). Với $m$ lớn hơn mà có ước nguyên tố lớn, vòng chia thử có thể tới cỡ $2·{10}^9$ bước; khi đó thay bước phân tích bằng Pollard rho.

Khi $m$ nguyên tố nhưng không lớn hơn $n$ (tức $p\le n$, ví dụ $p=13$), $n!$ chứa thừa số $p$ nên $n!\equiv 0$ và công thức giai thừa cũng hỏng: dùng định lý Lucas.

Khi nào phải đếm chính xác

Đề không cho modulo nghĩa là cần giá trị thật. Các mốc nên thuộc:

Vượt các mốc này thì dùng Pascal (chỉ cộng) trên số lớn, hoặc Python với số nguyên không giới hạn.

Ví dụ chạy tay

Tam giác Pascal tới hàng 5

 n\k   0   1   2   3   4   5
  0    1
  1    1   1
  2    1   2   1
  3    1   3   3   1
  4    1   4  [6]  4   1
  5    1   5  10  10   5   1
               ^
         C(5,2) = C(4,1) + C(4,2) = 4 + [6] = 10

Mỗi ô bằng ô ngay trên cộng ô trên bên trái. Lấy modulo sau mỗi phép cộng là đủ, với bất kỳ modulo nào.

Giai thừa và nghịch đảo giai thừa với $p=7$

Tính $\left(\genfrac{}{}{0ex}{}{5}{2}\right)mod7$ (đáp án thật $10\equiv 3$).

 i          :  0   1   2   3   4   5
 fact[i]    :  1   1   2   6   3   1      (24 = 3 mod 7, 120 = 1 mod 7)
 inv_fact[5] = fact[5]^(7-2) = 1^5 = 1
 đi lùi inv_fact[i-1] = inv_fact[i] * i :
 inv_fact   :  1   1   4   6   5   1
                       ^   ^
               kiểm tra: 2*4 = 8 = 1,  6*6 = 36 = 1  (mod 7)

 C(5,2) = fact[5] * inv_fact[2] * inv_fact[3] = 1 * 4 * 6 = 24 = 3 (mod 7)   đúng

Chia kẹo: 5 viên cho 3 em

 * * | * | * *      ->  (2, 1, 2)
 | * * * * * |      ->  (0, 5, 0)
 * * * * * | |      ->  (5, 0, 0)

 7 ô, chọn 2 ô làm vách:  C(7, 2) = 21 cách

Phản xạ với $n=3$ cặp ngoặc

Đường ())(() là đường xấu: nó xuống tới $-1$ ở bước 3. Lật phần từ đầu tới điểm chạm đầu tiên qua đường $y=-1$ (đổi ( thành ) và ngược lại, điểm xuất phát $0$ thành $-2$), giữ nguyên phần sau:

 bước       :   0   1   2   3   4   5   6
 đường xấu  :       (   )   )   (   (   )
 độ cao     :   0   1   0  -1   0   1   0
                           ^ chạm -1 lần đầu ở bước 3

 đường lật  :       )   (   (   (   (   )
 độ cao     :  -2  -3  -2  -1   0   1   0
                ^ xuất phát từ -2      phần sau bước 3 giữ nguyên

Ngược lại, mọi đường từ $-2$ về $0$ đều phải cắt $y=-1$, lật lại phần trước điểm cắt đầu tiên sẽ ra đúng một đường xấu. Đường sau khi lật đi từ $-2$ lên $0$ trong $6$ bước: $4$ bước lên, $2$ bước xuống. Số đường như vậy là $\left(\genfrac{}{}{0ex}{}{6}{4}\right)=15$, bằng số đường xấu. Vậy $C_3=\left(\genfrac{}{}{0ex}{}{6}{3}\right)-\left(\genfrac{}{}{0ex}{}{6}{4}\right)=20-15=5$, khớp năm dãy ((())), (()()), (())(), ()(()), ()()().

Legendre: $\left(\genfrac{}{}{0ex}{}{10}{4}\right)mod12$

$4!=24\equiv 0\left(\mathrm{mod}12\right)$ nên không nghịch đảo được. Tính số mũ từng nguyên tố:

 p  | v_p(10!)        | v_p(4!)   | v_p(6!)   | e_p
 ---+-----------------+-----------+-----------+-----
 2  | 5 + 2 + 1 = 8   | 2 + 1 = 3 | 3 + 1 = 4 |  1
 3  | 3 + 1     = 4   | 1         | 2         |  1
 5  | 2               | 0         | 1         |  1
 7  | 1               | 0         | 0         |  1

 C(10,4) = 2^1 * 3^1 * 5^1 * 7^1 = 210 = 6 (mod 12)

Cài đặt

Tam giác Pascal (modulo bất kỳ, $n$ nhỏ)

Chỉ cần một hàng nên dùng mảng một chiều: đi $j$ từ phải sang trái để C[j - 1] vẫn là giá trị của hàng trước. Nếu cần nhiều truy vấn $\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ khác nhau thì lưu cả bảng, nhưng nhớ bảng long long $2001\times 2001$ tốn khoảng $32$ MB (dùng int khi $M<2^{31}$ để giảm còn một nửa).

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

int main() {
    int n, k; long long M;
    cin >> n >> k >> M;                       // any modulus 1 <= M < 2^62
    if (k < 0 || k > n) { cout << 0 << "\n"; return 0; }
    vector<long long> C(n + 1, 0);            // one row of Pascal's triangle, O(n) memory
    C[0] = 1 % M;                             // row 0: C(0,0) = 1
    for (int i = 1; i <= n; i++)
        for (int j = i; j >= 1; j--)          // right to left: C[j-1] still holds row i-1
            C[j] = (C[j] + C[j - 1]) % M;     // Pascal: only additions, no division
    cout << C[k] << "\n";
}

Giai thừa + nghịch đảo giai thừa (modulo nguyên tố)

Khung chuẩn, dùng cho tổ hợp, chia kẹo và Catalan. Mỗi truy vấn $O(1)$.

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

const long long MOD = 1'000'000'007;          // prime
const int MAXN = 2'000'000;                   // must cover the largest n used
long long fact[MAXN + 1], inv_fact[MAXN + 1];

long long power(long long a, long long e) {   // a^e mod MOD, O(log e)
    long long r = 1; a %= MOD;
    while (e > 0) {
        if (e & 1) r = r * a % MOD;
        a = a * a % MOD;
        e >>= 1;
    }
    return r;
}

void precompute() {
    fact[0] = 1;
    for (int i = 1; i <= MAXN; i++) fact[i] = fact[i - 1] * i % MOD;
    inv_fact[MAXN] = power(fact[MAXN], MOD - 2);          // Fermat: a^(p-2) = a^(-1)
    for (int i = MAXN; i > 0; i--) inv_fact[i - 1] = inv_fact[i] * i % MOD;
}

long long C(long long n, long long k) {
    if (k < 0 || k > n) return 0;                           // out of range: 0 ways
    return fact[n] * inv_fact[k] % MOD * inv_fact[n - k] % MOD;
}

long long stars_bars(long long n, long long k) {            // n candies, k kids, >= 0 each
    if (k == 0) return n == 0 ? 1 : 0;                      // no kids: only n = 0 works
    return C(n + k - 1, k - 1);
}

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

// Input: q, then q lines, each one of
//   C n k  -> C(n, k)
//   S n k  -> n candies, k kids (>= 0 each)
//   K n    -> Catalan number C_n (only one number)
int main() {
    precompute();
    int q; cin >> q;
    while (q--) {
        char type; long long a, b = 0;
        cin >> type >> a;
        if (type != 'K') cin >> b;                              // 'K' has no second number
        if (type == 'C') cout << C(a, b) << "\n";               // binomial
        else if (type == 'S') cout << stars_bars(a, b) << "\n"; // stars and bars
        else cout << catalan(a) << "\n";                        // 'K': Catalan number
    }
}

Lưu ý MAXN phải phủ đối số lớn nhất thực sự dùng: chia kẹo cần tới $n+k-1$, Catalan cần tới $2n$. Hàm stars_bars tách riêng trường hợp $k=0$: công thức $\left(\genfrac{}{}{0ex}{}{n-1}{-1}\right)$ luôn trả $0$, trong khi theo quy ước chia $0$ viên cho $0$ em có đúng $1$ cách.

Dãy ngoặc khi biết trước tiền tố (phản xạ)

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

const long long MOD = 1'000'000'007;
vector<long long> fact, inv_fact;

long long power(long long a, long long e) {
    long long r = 1; a %= MOD;
    for (; e > 0; e >>= 1, a = a * a % MOD)
        if (e & 1) r = r * a % MOD;
    return r;
}
long long C(long long n, long long k) {
    if (k < 0 || k > n) return 0;
    return fact[n] * inv_fact[k] % MOD * inv_fact[n - k] % MOD;
}

int main() {
    int n; string s;
    cin >> n >> s;
    fact.assign(n + 1, 1); inv_fact.assign(n + 1, 1);
    for (int i = 1; i <= n; i++) fact[i] = fact[i - 1] * i % MOD;
    inv_fact[n] = power(fact[n], MOD - 2);
    for (int i = n; i > 0; i--) inv_fact[i - 1] = inv_fact[i] * i % MOD;

    int bal = 0;
    for (char c : s) {
        bal += (c == '(') ? 1 : -1;
        if (bal < 0) { cout << 0 << "\n"; return 0; }   // prefix already invalid
    }
    int r = n - (int)s.size();                          // characters still to place
    if ((r + bal) % 2 != 0 || bal > r) { cout << 0 << "\n"; return 0; }
    int d = (r + bal) / 2;                              // number of ')' still needed
    // all paths from height bal to 0, minus the bad ones (reflection)
    cout << (C(r, d) - C(r, d + 1) + MOD) % MOD << "\n";
}

$\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ theo modulo bất kỳ bằng công thức Legendre

Chương trình đầy đủ; main in $\left(\genfrac{}{}{0ex}{}{10}{4}\right)mod12=6$ như ví dụ chạy tay. Phép nhân đi qua __int128 nên đúng với mọi $1\le M<2^{62}$, kể cả $M=2^{32}$. Cách này không cần phân tích $M$, nên $M$ lớn cũng không làm chậm.

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

long long mul_mod(long long a, long long b, long long m) {   // (a*b) % m without overflow
    return (long long)((__int128)a * b % m);
}

long long power_mod(long long a, long long e, long long m) {
    long long r = 1 % m; a %= m;
    for (; e > 0; e >>= 1, a = mul_mod(a, a, m))
        if (e & 1) r = mul_mod(r, a, m);
    return r;
}

// exponent of prime p in x!  (Legendre)
long long legendre(long long x, long long p) {
    long long e = 0;
    while (x > 0) { x /= p; e += x; }
    return e;
}

// C(n,k) mod M for any 1 <= M < 2^62, n up to ~1e7
long long binom_any_mod(int n, int k, long long M) {
    if (k < 0 || k > n) return 0;
    vector<bool> composite(n + 1, false);            // sieve primes up to n
    long long ans = 1 % M;
    for (long long p = 2; p <= n; p++) {
        if (composite[p]) continue;
        for (long long q = p * p; q <= n; q += p) composite[q] = true;
        long long e = legendre(n, p) - legendre(k, p) - legendre(n - k, p);
        if (e > 0) ans = mul_mod(ans, power_mod(p, e, M), M);
    }
    return ans;
}

int main() {
    cout << binom_any_mod(10, 4, 12) << "\n";          // 6
    cout << binom_any_mod(17, 3, 1LL << 32) << "\n";   // 680
}

Cùng cách này tính được mọi biểu thức dạng tích và thương các giai thừa, ví dụ Catalan $C_m=\frac{(2m)!}{m!(m+1)!}$: chỉ cần đổi công thức $e_p$.

Nhiều truy vấn với modulo hợp số: tách ước nguyên tố của modulo

Chương trình đầy đủ, cũng nhân qua __int128 nên mọi phép nhân đúng với $1\le M<2^{62}$. Tuy vậy hàm dựng phân tích $M$ bằng chia thử tới $\sqrt{M}$, tốn $O(\sqrt{M})$: nhanh khi $M$ tới khoảng ${10}^{12}$ (cỡ ${10}^6$ bước), nhưng với $M$ gần $2^{62}$ có ước nguyên tố lớn thì vòng này tới cỡ $2·{10}^9$ bước, quá chậm. Khi đó thay đoạn chia thử bằng Pollard rho; phần còn lại giữ nguyên.

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

// Many queries C(n,k) mod M (M not necessarily prime).
// Multiplication is overflow-safe for M < 2^62, but the trial-division
// factoring below is fast only for M up to ~1e12 (use Pollard rho beyond).
// Every i = (part made of M's primes) * (part coprime to M).
// The coprime part is invertible mod M; the other part is kept as exponents.
struct BinomAnyMod {
    long long M;
    vector<long long> ps;                    // distinct primes of M
    vector<long long> fr;                    // fr[i] = product of coprime parts of 1..i, mod M
    vector<vector<int>> ex;                  // ex[j][i] = exponent of ps[j] in i!

    long long mul(long long a, long long b) const {   // (a*b) % M via __int128, no overflow
        return (long long)((__int128)a * b % M);
    }
    static long long ext_gcd(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, g = ext_gcd(b, a % b, x1, y1);
        x = y1; y = x1 - (a / b) * y1;
        return g;
    }
    long long inv(long long a) const {       // a must be coprime to M
        long long x, y;
        ext_gcd(a, M, x, y);
        return (x % M + M) % M;
    }
    long long power(long long a, long long e) const {
        long long r = 1 % M; a %= M;
        for (; e > 0; e >>= 1, a = mul(a, a))
            if (e & 1) r = mul(r, a);
        return r;
    }

    BinomAnyMod(int maxn, long long mod) : M(mod) {
        long long t = M;                     // factor M by trial division, O(sqrt M)
        for (long long p = 2; p * p <= t; p++)
            if (t % p == 0) { ps.push_back(p); while (t % p == 0) t /= p; }
        if (t > 1) ps.push_back(t);

        fr.assign(maxn + 1, 1 % M);
        ex.assign(ps.size(), vector<int>(maxn + 1, 0));
        for (int i = 1; i <= maxn; i++) {
            long long v = i;
            for (size_t j = 0; j < ps.size(); j++) {
                int c = 0;
                while (v % ps[j] == 0) { v /= ps[j]; c++; }
                ex[j][i] = ex[j][i - 1] + c;
            }
            fr[i] = mul(fr[i - 1], v % M);
        }
    }

    long long C(int n, int k) const {
        if (k < 0 || k > n) return 0;
        long long r = mul(mul(fr[n], inv(fr[k])), inv(fr[n - k]));
        for (size_t j = 0; j < ps.size(); j++)
            r = mul(r, power(ps[j], ex[j][n] - ex[j][k] - ex[j][n - k]));
        return r;
    }
};

int main() {
    BinomAnyMod small(100, 12);
    cout << small.C(10, 4) << "\n";          // 6, matches the hand-worked example
    BinomAnyMod big(100, 1LL << 32);
    cout << big.C(17, 3) << "\n";            // 680
}

Python: đếm chính xác và bản modulo

from math import comb

# exact values: Python ints never overflow
print(comb(100, 50))                 # 100891344545564193334812497256
print(comb(2 * 40, 40) // 41)        # Catalan C_40, already beyond long long

MOD = 10**9 + 7
N = 10**6
fact = [1] * (N + 1)
for i in range(1, N + 1):
    fact[i] = fact[i - 1] * i % MOD
inv_fact = [1] * (N + 1)
inv_fact[N] = pow(fact[N], MOD - 2, MOD)      # Fermat
for i in range(N, 0, -1):
    inv_fact[i - 1] = inv_fact[i] * i % MOD

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

math.comb tiện cho đếm chính xác với $n$ vài nghìn, nhưng không dùng nó rồi mới % MOD khi $n$ tới ${10}^6$: số trung gian có hàng trăm nghìn chữ số, rất chậm.

Độ phức tạp

⚠️ Lỗi thường gặp

• Tràn số khi nhân ba thừa số. fact[n] * inv_fact[k] * inv_fact[n-k] % MOD tràn long long vì tích hai số cỡ ${10}^9$ đã gần ${10}^{18}$. Phải % MOD sau mỗi phép nhân, như trong code mẫu.
• Modulo lớn hơn khoảng $3·{10}^9$. Khi $m^2\ge 2^{63}$ (ví dụ $m=2^{32}$), ngay cả a * b % m với $a,b<m$ cũng tràn long long. Nhân qua __int128: (__int128)a * b % m.
• Dùng Fermat khi modulo không nguyên tố. Với $m={10}^9$, $m=2^{32}$ hay $m$ đọc từ input, power(x, m - 2) cho kết quả sai mà không báo lỗi. Đọc kỹ modulo trong đề; nếu hợp số dùng Legendre, tách ước nguyên tố của modulo, hoặc Pascal.
• Modulo nguyên tố nhưng không lớn hơn $n$. Khi $p\le n$ thì fact[n] % p == 0, mọi $\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ tính ra đều $0$. Dùng Lucas.
• Mảng giai thừa không đủ dài. Chia kẹo gọi $\left(\genfrac{}{}{0ex}{}{n+k-1}{k-1}\right)$ và Catalan gọi $\left(\genfrac{}{}{0ex}{}{2n}{n}\right)$: kích thước cần là $n+k$ hay $2n$, không phải $n$. Thiếu thì truy cập ngoài mảng, kết quả rác.
• Hiệu hai tổ hợp bị âm. $\left(\genfrac{}{}{0ex}{}{2n}{n}\right)-\left(\genfrac{}{}{0ex}{}{2n}{n+1}\right)$ tính theo modulo có thể ra số âm; luôn (a - b + MOD) % MOD.
• Quên trường hợp $k<0$ hay $k>n$. Phản xạ có thể gọi C(r, d + 1) với $d+1>r$; chia kẹo với $x_i\ge 1$ có thể gọi C(n - 1, k - 1) với $n<k$. Hàm C phải trả $0$ ở các trường hợp này, nếu không sẽ đọc chỉ số âm.
• Chia cho số lần lặp không đều. "Đếm có thứ tự rồi chia" chỉ đúng khi mọi cấu hình bị đếm lặp như nhau. Ví dụ đếm vòng tròn xếp hạt có đối xứng (xoay, lật) thì số lần lặp khác nhau giữa các cấu hình: phải dùng Burnside.
• Cộng các trường hợp chồng lên nhau. "Có ít nhất một chữ số 7" không phải là cộng "vị trí 1 là 7" + "vị trí 2 là 7" + ...; các nhóm này giao nhau. Đếm phần bù (${10}^n-9^n$) hoặc bao hàm loại trừ.
• Đếm chính xác nhưng dùng long long. $\left(\genfrac{}{}{0ex}{}{67}{33}\right)$, $21!$, $C_{36}$ đã vượt $2^{63}$. Ngay cả khi đáp án vừa, công thức lặp r = r * (n - i) / (i + 1) vẫn có thể tràn ở bước nhân trung gian; dùng __int128 cho phép nhân hoặc Pascal.
• Mảng Pascal hai chiều quá lớn. long long C[2001][2001] đã khoảng $32$ MB, C[5001][5001] khoảng $200$ MB. Chỉ cần một hàng thì dùng mảng một chiều.

Biến thể / Mở rộng

• Định lý Lucas cho $\left(\genfrac{}{}{0ex}{}{n}{k}\right)modp$ khi $n$ tới ${10}^{18}$ và $p$ nhỏ; kết hợp CRT cho modulo hợp số không bình phương: Lucas Theorem.
• Bao hàm loại trừ: chia kẹo có chặn trên, đếm song ánh không điểm cố định (derangement), đếm theo tập con: Bao hàm loại trừ.
• Đếm dưới tác động đối xứng (vòng cổ, tô màu khối lập phương): Burnside và Pólya.
• Đếm có điều kiện phức tạp thường chuyển sang quy hoạch động đếm; công thức tổ hợp lúc đó đóng vai trò chuyển trạng thái nhanh, ví dụ QHĐ chữ số.
• Thứ hạng của hoán vị / tổ hợp (hoán vị thứ $k$, tổ hợp đứng thứ mấy theo từ điển) dùng giai thừa và $\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ để "nhảy" qua cả khối: Hạng hoán vị.
• Phân tích modulo lớn (tới cỡ ${10}^{18}$) để tách ước nguyên tố: Miller-Rabin và Pollard rho.
• Nền tảng số học đồng dư (lũy thừa nhanh, nghịch đảo, CRT): Toán học.

Bài tập luyện

Luyện kỹ thuật

• Hệ số nhị thức (binocoef) (Nghiệp dư) ${10}^5$ truy vấn $\left(\genfrac{}{}{0ex}{}{a}{b}\right)mod{10}^9+7$ với $a\le {10}^6$: bài luyện trực tiếp khung giai thừa + nghịch đảo giai thừa.
• Phân phát táo (distapples) (Nghiệp dư) chia $m$ quả táo giống nhau cho $n$ em, mỗi em được $0$ quả trở lên: đúng bài chia kẹo, chú ý kích thước mảng giai thừa.
• Dãy ngoặc I (bracketseq1) (Nghiệp dư) đếm dãy ngoặc đúng độ dài $n$: số Catalan, nhớ xử lý $n$ lẻ.
• Dãy ngoặc II (bracketseq2) (Trung cấp) đã biết một tiền tố: áp dụng nguyên lý phản xạ từ độ cao hiện tại thay vì từ $0$.
• Dãy điểm của Dima (twoseq) (Nâng cao) đáp án là tích các giai thừa chia cho lũy thừa của $2$, lấy theo modulo $m$ tùy ý (có thể chẵn): không nghịch đảo được, phải xử lý thừa số $2$ riêng.

Đề VOI

• Kế hoạch cải tổ (voi15_reform) (Nâng cao) VOI 2015: phân loại cạnh bị bỏ theo cầu / không cầu, rồi đếm từng loại bằng quy tắc cộng và quy tắc nhân.
• Phép toán OR (voi21_or) (Nâng cao) VOI 2021: số bộ $k$ chỉ số chọn trong $c$ chỉ số "được phép" là $\left(\genfrac{}{}{0ex}{}{c}{k}\right)mod{10}^9+7$; cần bảng giai thừa, sau đó kết hợp với bao hàm loại trừ trên tập con.
• Hàng cây (voi11_tline) (Kỳ cựu) VOI 2011: câu hai đếm hoán vị không có 3 phần tử giảm dần; liệt kê đáp án với số cây nhỏ để nhận ra dãy quen thuộc, và để ý modulo ${10}^9$ không nguyên tố.
• Sản xuất đồ chơi (voi13_organ) (Kỳ cựu) VOI 2013: số cách xếp một lô chỉ phụ thuộc số ống, lập bảng truy hồi kiểu Pascal bằng cách chèn ống cao nhất; đề không lấy modulo nên cần số lớn.
• Tất niên (voi26_tet) (Kỳ cựu) VOI 2026: độ đa dạng của một đoạn: đếm các cặp vị trí có giá trị khác nhau, mỗi cặp cho đúng một dãy mới; phần khó còn lại là duy trì đại lượng này khi đoạn thay đổi.
• Mã hóa dãy số (voi25_encode) (Chuyên gia) VOI 2025: sau khi biết có bao nhiêu vị trí buộc phải tăng ngặt, bài toán còn lại là đếm dãy không giảm với giá trị bị chặn, tức chia kẹo theo modulo ${10}^9+7$.