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 phần tử khác nhau (không phân biệt thứ tự) của tập hợp phần tử đã cho () được gọi là một tổ hợp chập của phần tử.
Số các tổ hợp chập của phần tử khác nhau được kí hiệu là hoặc :
Để 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 .
Ta quy ước:
Với công thức này, ta nghĩ ngay đến một thuật toán "ngây thơ": Tính , và . Từ đó tính được .
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:
Vì là số nguyên, nên bạn yên tâm rằng luôn chia hết cho .
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 có thể rất lớn (khi thì )
Vì kết quả của 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 theo một modulo nào đó (phần dư của phép chia cho ).
Dưới đây là một số cách sử dụng để tính theo modulo . Chú ý rằng độ phức tạp dưới đây được tính toán khi sử dụng một modulo 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 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 đó:
Cho là một số nguyên tố và số nguyên không chia hết cho . Khi đó, ta có:
Từ đó, ta rút ra:
Lũy thừa nhanh
Ý tưởng:
Ta viết lại:
Ta sử dụng hai mảng: mảng để lưu và mảng để lưu . Sau đó dùng công thức (sử dụng định lý Fermat nhỏ):
Chú ý rằng nên ta chỉ tính và với .
Ta sẽ tính mảng như sau:
ế
Tiếp theo ta sử dụng thuật toán lũy thừa nhanh để tính với độ phức tạp .
Còn mảng thì tính như sau:
ế
Cuối cùng, .
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';
}
}
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';
}
}
Định lý Euler (mở rộng của định lý Fermat nhỏ)
Cho số nguyên nguyên tố cùng nhau. Khi đó, ta có: .
Trong đó, là hàm phi Euler: êố.
Từ đó, ta rút ra:
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)
Trong đó:
ếàòạ
Bạn đọc có thể thấy, là số mũ của khi phân tích ra thừa số nguyên tố.
là tích tất cả các số từ đến và không bao gồm các số chia hết cho (với là số nguyên tố).
định nghĩa tương tự
là vị trí cuối cùng mà . Nghĩa là ta chỉ cần chạy cho đến khi .
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ị, sử dụng để lưu (Xem phần mô tả mở rộng định lý Lucas) thay vì 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;
}
Định lý thặng dư Trung Hoa là cầu nối giúp ta tính toán được khi 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ệ:
với đôi một nguyên tố cùng nhau.
Ký hiệu:
Từ đó nhận thấy:
Khi này chỉ cần cộng tất cả số hạng ta được một nghiệm thỏa mãn hệ:
Để 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 với
Lời giải:
Giả sử bằng những cách trên, bạn đã tính được 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
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 trong công thức để 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();
}