Người viết:
- Nguyễn Minh Hiển – Trường Đại học Công nghệ, ĐHQGHN
Reviewer:
- Phạm Công Minh – Trường THPT chuyên Khoa học Tự nhiên, ĐHQGHN
- Nguyễn Đức Kiên – Trường Đại học Công nghệ, ĐHQGHN
- Nguyễn Minh Nhật – Trường THPT chuyên Khoa học Tự nhiên, ĐHQGHN
¶ Tổ hợp (Combination)
¶ Giới thiệu về tổ hợp
Trong toán học, tổ hợp là cách chọn các phần tử từ một nhóm mà không phân biệt thứ tự chọn. Mỗi tập con gồm kkk phần tử khác nhau (không phân biệt thứ tự) của tập hợp nnn phần tử đã cho (0≤k≤n0 ≤ k ≤ n0≤k≤n) được gọi là một tổ hợp chập kkk của nnn phần tử. Số các tổ hợp chập kkk của nnn phần tử khác nhau được kí hiệu là CnkC_n^kCnk hoặc (nk)dbinom{n}{k}(kn):
Cnk=(nk)=n!k!(n−k)!=n(n−1)…(n−k+1)k!, 0≤k≤nC_n^k = dbinom{n}{k} = dfrac{n!}{k! (n – k)!} = dfrac{n(n-1)ldots(n – k + 1)}{k!},, 0 ≤ k ≤ n Cnk=(kn)=k!(n−k)!n!=k!n(n−1)…(n−k+1),0≤k≤n
Để bạn đọc tiện theo dõi, trong bài viết này, chúng ta thống nhất sử dụng ký hiệu CnkC_n^kCnk. Ta quy ước:
- Cn0=1C_n^0 = 1Cn0=1 với mọi n≥0n ge 0n≥0
- Nếu k>nk > nk>n thì Cnk=0C_n^k = 0Cnk=0
¶ Một số tính chất của tổ hợp
-
Với mọi n≥1n ≥ 1n≥1 và 0≤k≤n0 ≤ k ≤ n0≤k≤n, ta có:
- Cnk=Cnn−kC_n^k = C_n^{n – k}Cnk=Cnn−k
- Cnk=Cn−1k−1+Cn−1kC_n^k = C_{n – 1}^{k – 1} + C_{n – 1}^{k}Cnk=Cn−1k−1+Cn−1k
-
CnkC_n^kCnk còn được gọi là hệ số nhị thức (binomial coefficients) do CnkC_n^kCnk là hệ số trong khai triển:
(x+y)n=∑k=0nCnk⋅xk⋅yn−k(x + y)^n= sumlimits_{k = 0}^{n} C_n^k cdot x^k cdot y ^ {n – k} (x+y)n=k=0∑nCnk⋅xk⋅yn−k
¶ Tính số tổ hợp
¶ Sử dụng định nghĩa
Cnk=n!k!(n−k)!C_n^k = dfrac{n!}{k! (n – k)!} Cnk=k!(n−k)!n!
- Với công thức này, ta nghĩ ngay đến một thuật toán “ngây thơ”: Tính n!n!n!, k!k!k! và (n−k)!(n – k)!(n−k)!. Từ đó tính được CnkC_n^kCnk.
long long res = 1; for (int i = 1; i <= n; i++) res = res * i; for (int i = 1; i <= k; i++) res = res / i; for (int i = 1; i <= n-k; i++) res = res / i;
- Mở rộng hơn, ta có thể biến đổi một chút như sau:
Cnk=n1⋅n−12…n−k+1k=Cn−1k−1⋅(n−k+1)kC_n^k = dfrac{n}{1} cdot dfrac{n – 1}{2} ldots dfrac{n – k + 1}{k} = dfrac{C_{n-1}^{k-1} cdot (n – k + 1)}{k} Cnk=1n⋅2n−1…kn−k+1=kCn−1k−1⋅(n−k+1)
Vì CnkC_n^kCnk là số nguyên, nên bạn yên tâm rằng Cn−1k−1⋅(n−k+1)C_{n-1}^{k-1} cdot (n – k + 1)Cn−1k−1⋅(n−k+1) luôn chia hết cho kkk.
long long res = 1; for (int i = 1; i <= k; i++) res = res * (n – i + 1) / i;
- Hai cách tiếp cận trên rất tự nhiên, dễ nghĩ, dễ thực hiện nhưng lại có một trở ngại: giá trị của n!n!n! có thể rất lớn (khi n=20n = 20n=20 thì n!≈2.42×1018n! approx 2.42times10^{18}n!≈2.42×1018)
¶ Sử dụng công thức truy hồi
Cnk={1neˆˊu k=0 vaˋ n≥00neˆˊu k>nCn−1k−1+Cn−1ktrong caˊc trường hợp coˋn lạibegin{align*} C_n^k=begin{cases} 1 &text{nếu } k = 0 text{ và } n ge 0 0 &text{nếu } k > n C_{n – 1}^{k – 1} + C_{n – 1}^{k} &text{trong các trường hợp còn lại} end{cases} end{align*}Cnk=⎩⎨⎧10Cn−1k−1+Cn−1kneˆˊu k=0 vaˋ n≥0neˆˊu k>ntrong caˊc trường hợp coˋn lại
Với công thức truy hồi này, ta sẽ sử dụng một mảng hai chiều C[n][k] để tính CnkC_n^kCnk
Code C++ minh họa
for (int i = 0; i <= n; i++){ C[i][0] = 1; for (int k = 1; k <= i; k++){ C[i][k] = C[i – 1][k – 1] + C[i – 1][k]; } }
Độ phức tạp không gian: O(n2)O(n^2)O(n2) Độ phức tạp thời gian:
- Tiền xử lý: O(n2)O(n^2)O(n2)
- Truy vấn: O(1)O(1)O(1)
¶ Tính số tổ hợp theo modulo M
Vì kết quả của CnkC_n^kCnk có thể rất lớn nên trong lập trình thi đấu, thí sinh thường được yêu cầu tính CnkC_n^kCnk theo một modulo MMM nào đó (phần dư của phép chia CnkC_n^kCnk cho MMM). Dưới đây là một số cách sử dụng để tính CnkC_n^kCnk theo modulo MMM. Chú ý rằng độ phức tạp dưới đây được tính toán khi sử dụng một modulo MMM cho toàn bài!
Ngoài ra còn có hai cách tính dựa trên cách tính giai thừa modulo MMM khá hiệu quả. Tham khảo thêm tại Giai thừa modulo p . Dưới đây là đánh giá về hai cách đó:
¶ Sử dụng công thức truy hồi
Ở đây, ta sẽ sử dụng công thức truy hồi ở trên và thay đổi một chút:
Cnk=(Cn−1k−1+Cn−1k)mod MC_n^k = (C_{n – 1}^{k – 1} + C_{n – 1}^{k}) mod M Cnk=(Cn−1k−1+Cn−1k)modM
Code C++ minh họa
for (int i = 0; i <= n; i++){ C[i][0] = 1 % MOD; for (int k = 1; k <= i; k++){ C[i][k] = (C[i – 1][k – 1] + C[i – 1][k]) % MOD; } }
Độ phức tạp không gian: O(n2)O(n^2)O(n2) Độ phức tạp thời gian:
- Tiền xử lý: O(n2)O(n^2)O(n2)
- Truy vấn: O(1)O(1)O(1)
Nhận xét: Đây là cách đơn giản, dễ nghĩ, dễ code đúng, nên sử dụng trong trường hợp nnn nhỏ để tiết kiệm thời gian.
¶ Sử dụng định nghĩa
Rào cản lớn nhất cho việc sử dụng định nghĩa Cnk=n!k!(n−k)!C_n^k = dfrac{n!}{k! (n – k)!}Cnk=k!(n−k)!n! là n!n!n! quá lớn. Tuy nhiên khi ta cần lấy kết quả theo modulo MMM, đó lại là vấn đề khác.
-
Điều kiện sử dụng: MMM nguyên tố và n<Mn < Mn<M
-
Kiến thức sử dụng:
- Nghịch đảo modulo (Modular Inverse): Tham khảo ở bài viết Số học 4.5 – Nghịch đảo modulo của VNOI.
- Định lý Fermat nhỏ
Chứng minh
Cho ppp là một số nguyên tố và số nguyên aaa không chia hết cho ppp. Khi đó, ta có:
ap−1≡1(modp)a^{p – 1} equiv 1 pmod p ap−1≡1(modp)
Từ đó, ta rút ra:
a−1≡ap−2(modp)a^{-1} equiv a^{p-2} pmod p a−1≡ap−2(modp)
- Lũy thừa nhanh
-
Ý tưởng:
- Ta viết lại:
Cnk=n!×(k!)−1×((n−k)!)−1mod MC_n^k = n! times left( k! right)^{-1} times left( (n – k)! right)^{-1} mod M Cnk=n!×(k!)−1×((n−k)!)−1modM
-
Ta sử dụng hai mảng: mảng fact[i]text{fact}[i]fact[i] để lưu i! mod Mi! bmod Mi!modM và mảng ifact[i]text{ifact}[i]ifact[i] để lưu (i!)−1 mod M(i!)^{-1} bmod M(i!)−1modM. Sau đó dùng công thức (sử dụng định lý Fermat nhỏ):
ifact[i]=(fact[i])−1 mod M=(fact[i])M−2 mod Mtext{ifact}[i] = (text{fact}[i]) ^ {-1} bmod M = (text{fact}[i])^{M-2} bmod M ifact[i]=(fact[i])−1modM=(fact[i])M−2modM
Chú ý rằng fact[i]≡0(modM) ∀i≥Mtext{fact}[i] equiv 0 pmod M ;;forall i ge Mfact[i]≡0(modM)∀i≥M nên ta chỉ tính fact[i]text{fact}[i]fact[i] và ifact[i]text{ifact}[i]ifact[i] với 0≤i≤M−10 le i le M – 10≤i≤M−1.
-
Ta sẽ tính mảng fact[i]text{fact}[i]fact[i] như sau:
{fact[0]=1fact[i]=(fact[i−1]×i) mod M neˆˊu 1≤i≤nbegin{align} begin{cases} text{fact}[0] &= 1 text{fact}[i] &= (text{fact}[i – 1] times i ) bmod M &text{ nếu } 1 le i le n end{cases} end{align}{fact[0]fact[i]=1=(fact[i−1]×i)modM neˆˊu 1≤i≤n
- Tiếp theo ta sử dụng thuật toán lũy thừa nhanh để tính ifact[n]=(fact[n])M−2mod Mtext{ifact}[n] = left( text{fact}[n] right)^{M-2} mod Mifact[n]=(fact[n])M−2modM với độ phức tạp O(logM)O(log M)O(logM).
- Còn mảng ifact[i]text{ifact}[i]ifact[i] thì tính như sau:
{ifact[n]=(fact[n])M−2mod Mifact[i−1]=ifact[i]×imod Mneˆˊu 1≤i≤nbegin{align} begin{cases} text{ifact}[n] &= left( text{fact}[n] right)^{M-2} &mod M text{ifact}[i – 1] &= text{ifact}[i] times i &mod M &text{nếu } 1 le i le n end{cases} end{align}{ifact[n]ifact[i−1]=(fact[n])M−2=ifact[i]×imodMmodMneˆˊu 1≤i≤n
- Cuối cùng, Cnk=fact[n]×ifact[k]×ifact[n−k]mod MC_n^k = text{fact}[n] times text{ifact}[k] times text{ifact}[n – k] mod MCnk=fact[n]×ifact[k]×ifact[n−k]modM.
Code C++ minh họa
const int MOD = 1e9 + 7; const int N = 1e6; int fact[N + 5], ifact[N + 5]; // Hàm lũy thừa nhanh long long binpow(long long a, long long b) { long long ans = 1; while (b > 0){ if (b % 2) ans = ans * a % MOD; a = a * a % MOD; b /= 2; } return ans; } // Chuẩn bị void prepare(){ // Tính fact[] fact[0] = 1; for (int i = 1; i <= N; i++) fact[i] = 1LL * fact[i – 1] * i % MOD; // Tính ifact[] ifact[N] = binpow(fact[N], MOD – 2); for (int i = N – 1; i >= 1; i-) ifact[i] = 1LL * ifact[i + 1] * (i + 1) % MOD; } // Hàm tính nCk int C(int n, int k){ if (k > n) return 0; return (1LL * fact[n] * ifact[k] % MOD) * ifact[n – k] % MOD; } int main(){ prepare(); // Truy vấn int q; cin >> q; while (q-){ int n, k; cin >> n >> k; cout << C(n, k) << ‘n’; } }
Độ phức tạp không gian: O(n)O(n)O(n) Độ phức tạp thời gian:
- Tiền xử lý: O(n+logM)O(n + log M)O(n+logM)
- Truy vấn: O(1)O(1)O(1)
¶ Sử dụng định lý Lucas
Đây là cải tiến của cách sử dụng định nghĩa để có thể sử dụng cho n>Mn > Mn>M
- Điều kiện sử dụng: MMM nguyên tố
- Ý tưởng: Các bạn có thể tham khảo tại bài viết Định lý Lucas của VNOI
Code C++ minh họa
int C(long long n, long long k){…} // hàm tính Ckn sử dụng định nghĩa bên trên int comb(long long n, long long k){ if (k > n) return 0; int res = 1; while (n > 0){ res = 1LL * res * C(n % MOD, k % MOD) % MOD; n /= MOD; k/= MOD; } return res; } Bạn đọc tham khảo thêm code đầy đủ ở đây const int MOD = 1e6 + 3; int fact[MOD + 5], ifact[MOD + 5]; // Hàm lũy thừa nhanh long long binpow(long long a, long long b) { long long ans = 1; while (b > 0){ if (b % 2) ans = ans * a % MOD; a = a * a % MOD; b /= 2; } return ans; } // Chuẩn bị void prepare(){ // Tính fact[] fact[0] = 1; for (int i = 1; i < MOD; i++) fact[i] = 1LL * fact[i – 1] * i % MOD; // Tính ifact[] ifact[MOD – 1] = binpow(fact[MOD – 1], MOD – 2); for (int i = MOD – 2; i >= 0; i-) ifact[i] = 1LL * ifact[i + 1] * (i + 1) % MOD; } // Hàm tính nCk với n < M int C(int n, int k){ if (k > n) return 0; return (1LL * fact[n] * ifact[k] % MOD) * ifact[n – k] % MOD; } // Hàm tính nCk với n có thể lớn hơn M int comb(long long n, long long k){ if (k > n) return 0; int res = 1; while (n > 0){ res = 1LL * res * C(n % MOD, k % MOD) % MOD; n /= MOD; k/= MOD; } return res; } int main(){ prepare(); // Truy vấn int q; cin >> q; while (q-){ long long n, k; cin >> n >> k; cout << comb(n, k) << ‘n’; } }
Độ phức tạp không gian: O(M)O(M)O(M) Độ phức tạp thời gian:
- Tiền xử lý: O(M)O(M)O(M)
- Truy vấn: O(logM(n))Obig(log_M(n)big)O(logM(n))
¶ Sử dụng định lý Lucas mở rộng
- Điều kiện sử dụng: MMM lũy thừa của số nguyên tố
- Kiến thức sử dụng:
-
Nghịch đảo modulo, lũy thừa nhanh
-
Định lý Euler (mở rộng của định lý Fermat nhỏ) Cho 222 số nguyên a,ma, ma,m nguyên tố cùng nhau. Khi đó, ta có: aφ(m)≡1(modm)a^{varphi(m)} equiv 1 pmod maφ(m)≡1(modm).
Trong đó, φ(m)varphi(m)φ(m) là hàm phi Euler: φ(m)=m⋅∏p nguyeˆn toˆˊ,p∣mp−1pvarphi(m) = m cdot prodlimits_{p text{ nguyên tố}, p mid m} dfrac{p – 1}{p}φ(m)=m⋅p nguyeˆn toˆˊ,p∣m∏pp−1. Từ đó, ta rút ra:
a−1≡aφ(m)−1(modp)a^{-1} equiv a^{varphi(m)-1} pmod p a−1≡aφ(m)−1(modp)
-
Mở rộng định lý Lucas cho modulo là lũy thừa số nguyên tố: Andrew Granville đã chứng minh được công thức sau: (Xem bài báo tại đây hoặc tại đây)
teqpe1Cnk≡(n0!)p(k0!)p(r0!)p⋅(n1!)p(k1!)p(r1!)p…(nd!)p(kd!)p(rd!)p ( mod pq)dfrac{t^{e_q}}{p^{e_1}} C_n^k equiv dfrac{(n_0!)_p}{(k_0!)_p(r_0!)_p} cdot dfrac{(n_1!)_p}{(k_1!)_p(r_1!)_p} ldots dfrac{(n_d!)_p}{(k_d!)_p(r_d!)_p} ,(bmod p^q) pe1teqCnk≡(k0!)p(r0!)p(n0!)p⋅(k1!)p(r1!)p(n1!)p…(kd!)p(rd!)p(nd!)p(modpq)
Trong đó:
- t = begin{align} begin{cases} 1 &text{ nếu } p = 2 text{ và } q ge 3 -1&text{ còn lại} end{cases} end{align}
- ej=∑i≥j(⌊npi⌋−⌊kpi⌋−⌊rpi⌋)e_j = sumlimits_{i ge j} left( leftlfloor dfrac{n}{p^i} rightrfloor – leftlfloor dfrac{k}{p^i} rightrfloor – leftlfloor dfrac{r}{p^i} rightrfloor right)ej=i≥j∑(⌊pin⌋−⌊pik⌋−⌊pir⌋) Bạn đọc có thể thấy, e1e_1e1 là số mũ của ppp khi phân tích CnkC_n^kCnk ra thừa số nguyên tố.
- (n!)pleft( n! right)_p(n!)p là tích tất cả các số từ 111 đến nnn và không bao gồm các số chia hết cho ppp (với ppp là số nguyên tố).
- r=n−kr = n-kr=n−k
- ni=⌊npi⌋ mod pqn_i = leftlfloor dfrac{n}{p^i} rightrfloor bmod p ^ qni=⌊pin⌋modpq
- ki,rik_i, r_iki,ri định nghĩa tương tự nin_ini
- ddd là vị trí cuối cùng mà ni≠0n_i neq 0ni=0. Nghĩa là ta chỉ cần chạy cho đến khi ni=0n_i = 0ni=0.
-
Chi tiết các bước:
const int MOD = 27; const int prime = 3; long long fact[MOD], ifact[MOD];
- Chú ý rằng ở bước chuẩn bị, fact[i]text{fact}[i]fact[i] sử dụng để lưu (i!)pleft( i! right)_p(i!)p (Xem phần mô tả mở rộng định lý Lucas) thay vì i!i!i! và ta cần sử dụng định lý Euler thay cho định lý Fermat nhỏ.
void init(){ fact[0] = 1; for (int i = 1; i < MOD; i++) { if (i % prime != 0) fact[i] = (fact[i – 1] * i) % MOD; else fact[i] = fact[i – 1]; } int phi = MOD / prime * (prime – 1) – 1; ifact[MOD – 1] = binpow(fact[MOD – 1], phi, MOD); for (int i = MOD – 1; i > 0; i-) { if (i % prime != 0) ifact[i – 1] = (ifact[i] * i) % MOD; else ifact[i – 1] = ifact[i]; } }
- Tiếp theo ta sử dụng công thức bên trên
long long C(long long N, long long K, long long R){ return (fact[N] * ifact[R] % MOD) * ifact[K] % MOD; } int count_carry(long long n, long long k, long long r, int p, long long t){ long long res = 0; while (n >= t) { res += ((n / t) – (k / t) – (r / t)); t *= p; } return res; } long long calc(long long N, long long K, long long R) { if (K > N) return 0; long long res = 1; int vp = count_carry(N, K, R, prime, prime); int vp2 = count_carry(N, K, R, prime, MOD); while (N > 0) { res = (res * C(N % MOD, K % MOD, R % MOD)) % MOD; N /= prime; K /= prime; R /= prime; } res = res * binpow(prime, vp, MOD) % MOD; if ((vp2 % 2 == 1) && (prime != 2 || MOD <= 4)) res = (MOD – res) % MOD; return res; }
Độ phức tạp không gian: O(M)O(M)O(M) Độ phức tạp thời gian:
- Tiền xử lý: O(M)O(M)O(M)
- Truy vấn: O(logp(n))Obig(log_p(n)big)O(logp(n))
¶ Sử dụng định lý thặng dư Trung Hoa
Định lý thặng dư Trung Hoa là cầu nối giúp ta tính toán được khi MMM không phải là số nguyên tố.
- Kiến thức sử dụng:
-
Định lý Thặng dư trung hoa – Chinese remainder theorem (CRT)
Xét hệ:
{a≡a1(modm1)a≡a2(modm2)⋮a≡ak(modmk)left{begin{array}{rcl} a & equiv & a_1 pmod{m_1} a & equiv & a_2 pmod{m_2} & vdots & a & equiv & a_k pmod{m_k} end{array}right.⎩⎨⎧aaa≡≡⋮≡a1(modm1)a2(modm2)ak(modmk)
với m1,m2,…mkm_1, m_2, ldots m_km1,m2,…mk đôi một nguyên tố cùng nhau.
Ký hiệu:
-
M=m1⋅m2…mkM = m_1 cdot m_2 ldots m_kM=m1⋅m2…mk
-
Mi=MmiM_i = dfrac{M}{m_i}Mi=miM
-
Ni=Mi−1mod miN_i = M_i^{-1} mod m_iNi=Mi−1modmi
Từ đó nhận thấy:
{aiMiNi≡ai(modmi)ajMjNj≡0(modmi) ∀j≠ileft{begin{array}{rcl} a_i M_i N_i & equiv & a_i pmod{m_i} a_j M_j N_j & equiv & 0 pmod{m_i} ,,forall j neq i end{array}right.{aiMiNiajMjNj≡≡ai(modmi)0(modmi)∀j=i
Khi này chỉ cần cộng tất cả số hạng aiMiNia_i M_i N_iaiMiNi ta được một nghiệm thỏa mãn hệ:
a=a1M1N1+a2M2N2+…+akMkNka = a_1 M_1 N_1 + a_2 M_2 N_2 + ldots + a_k M_k N_k a=a1M1N1+a2M2N2+…+akMkNk
-
-
Để minh họa rõ hơn này, ta sẽ giải quyết bài nCr trên Hackerrank
Tóm tắt: Tính Cnk mod 142857C_n^k bmod 142857Cnkmod142857 với 0≤k≤n≤1090 le k le n le 10^{9}0≤k≤n≤109
Lời giải:
-
Giả sử bằng những cách trên, bạn đã tính được CnkC_n^kCnk modulo là số nguyên tố (ĐL Lucas) hoặc lũy thừa của chúng (ĐL Lucas mở rộng). Tiếp theo ta sẽ sử dụng CRT xử lý các phần còn lại.
-
Đầu tiên, ta sẽ phân tích modulo 142857=33⋅11⋅13⋅17142857 = 3^3 cdot 11 cdot 13 cdot 17142857=33⋅11⋅13⋅17
int n_primes = 4; int primes[] = {3, 11, 13, 37}; int primes_pw[] = {27, 11, 13, 37}; int rem[n_primes]; vector<long long> fact[n_primes], ifact[n_primes];
-
Ta chuẩn bị sẵn một mảng tính MiNiM_i N_iMiNi trong công thức a=∑aiMiNia = sum a_i M_i N_ia=∑aiMiNi để tiện cho việc truy vấn.
void prepare(){ for (int i = 0; i < n_primes; i++) { // M_i int tmp = MOD / primes_pw[i]; //giá trị phi Euler của primes_pw[i] int phi = primes_pw[i] / primes[i] * (primes[i] – 1); // M_i * M_i^(-1) rem[i] = tmp * binpow(tmp, phi – 1, primes_pw[i]) % MOD; } }
-
Cuối cùng tính ra kết quả cuối cùng sử dụng CRT
long long CRT(long long N, long long K) { long long res = 0; for (int i = 0; i < n_primes; i++) { // C(n, k, MOD) là hàm tính nCk modulo MOD int ans = C(N, K, MOD) * rem[i] % MOD; res = (res + ans) % MOD; } return res; }
Bạn đọc tham khảo code nộp AC bài nCr ở đây #include <bits/stdc++.h> const int MOD = 142857; using ll = long long; using namespace std; int primes[] = {3, 11, 13, 37}; int primes_pw[] = {27, 11, 13, 37}; int phi[] = {18, 10, 12, 36}; // phi = prime_pw * (prime – 1)/prime int rem[4]; vector<ll> fact[4], ifact[4]; int t; ll binpow(ll a, ll n, ll mod) { ll res = 1; for (; n > 0; n >>= 1) { if (n & 1) res = res * a % mod; a = a * a % mod; } return res; } void init(int x) { fact[x].assign(primes_pw[x], 0); ifact[x].assign(primes_pw[x], 0); fact[x][0] = 1; for (int i = 1; i < primes_pw[x]; i++) { if (i % primes[x] != 0) fact[x][i] = (fact[x][i – 1] * i) % primes_pw[x]; else fact[x][i] = fact[x][i – 1]; } ifact[x][primes_pw[x] – 1] = binpow(fact[x][primes_pw[x] – 1], primes_pw[x] / primes[x] * (primes[x] – 1) – 1, primes_pw[x]); for (int i = primes_pw[x] – 1; i > 0; i-) { if (i % primes[x] != 0) ifact[x][i – 1] = (ifact[x][i] * i) % primes_pw[x]; else ifact[x][i – 1] = ifact[x][i]; } } /*i is the order of prime*/ ll C(ll N, ll K, ll R, int i) { return (fact[i][N] * ifact[i][R] % primes_pw[i]) * ifact[i][K] % primes_pw[i]; } int count_carry(ll n, ll k, ll r, ll p, ll t) { ll res = 0; while (n >= t) { res += (n / t – k / t – r / t); t *= p; } return res; } ll calc(ll N, ll K, ll R, int ord_pr) { if (K > N) return 0; int prime = primes[ord_pr]; int mod = primes_pw[ord_pr]; ll res = 1; int vp = count_carry(N, K, R, prime, prime); int vp2 = count_carry(N, K, R, prime, mod); while (N > 0) { res = (res * C(N % mod, K % mod, R % mod, ord_pr)) % mod; N /= prime; K /= prime; R /= prime; } res = res * binpow(prime, vp, mod) % mod; if ((vp2 & 1) && (prime != 2 || mod <= 4)) res = (mod – res) % mod; return res; } ll CRT(ll N, ll K) { ll res = 0; for (int i = 0; i <= 3; i++) { int ans = calc(N, K, N – K, i) * rem[i] % MOD; res = (res + ans) % MOD; } return res; } void solve() { for (int i = 0; i <= 3; i++) { init(i); int tmp = MOD / primes_pw[i]; rem[i] = tmp * binpow(tmp, phi[i] – 1, primes_pw[i]) % MOD; } while (t-) { ll N, K; cin >> N >> K; cout << CRT(N, K) << ‘n’; } } int main() { ios_base::sync_with_stdio(0); cin.tie(NULL); cin >> t; solve(); }
¶ Bài tập luyện tập
- SPOJ – Marbles
- Codeforces 181C Div.2
- Codechef – Long Sandwich
- Hackerearth – Binomial Coefficient
- SPOJ – ADATEAMS
- SPOJ – UCV2013E
- Codeforces 239C Div. 1
- VNOJ – Xông nhà ngày Tết
- LQDOJ – Tổ hợp nCk
- LQDOJ – Tổ hợp nCk 1
- LQDOJ – Tổ hợp nCk 2
- LQDOJ – Tổ hợp nCk 3
- LQDOJ – Hệ số nhị thức
- Hackerrank – nCr
- Hackerearth – Army parade
- Library Checker – Binomial Coefficient
- Libary Checker – Binomial Coefficient (Prime Mod)
¶ Tài liệu tham khảo
- Codeforces – Binomial Coefficients with Large Mod
- Codeforces – nCk for small k
- Mở rộng định lý Lucas của Andrew Granville
- CP-Algorithms – Binomial Coefficients