如何求組合數\(C_a^b\)
一、預處理法一
例題：https://www.acwing.com/problem/content/887/
理論依據：\(\huge C_a^b=C_{a-1}^b+C_{a-1}^{b-1}\)
適合場景：
1、\(\large a<=2000,b<=2000\)
2、詢問次數：\(\large n<=1e5\)
我們有\(a\)個蘋果，現在需要選出\(b\)個蘋果。一共有多少種選法呢？
我們走到第一個蘋果面前，我們面臨兩個選擇：
（1)選擇它 （2）放棄它
如果選擇了它，將要面對\(a-1\)個蘋果中選\(b-1\)個蘋果的問題。
如果放棄它，將要面對\(a-1\)個蘋果中選\(b\)個蘋果的問題。
根據加法原理，兩種選擇方法加在一起，就是方法總數，就是這面的遞推式。
#include <bits/stdc++.h>
using namespace std;
const int N = 2010; //題目數據上限是2000，這里加了10
const int MOD = 1e9 + 7; //取模的質數
int c[N][N]; //結果數組
//排列組合 C(A,B)=C(A-1,B)+C(A-1,B-1)
void init() {
for (int i = 0; i < N; i++)//每個數都來計算
for (int j = 0; j <= i; j++)//j需要比i小~
//特判，從i個蘋果里選擇0個蘋果，只有一種方案，就是，啥也不選。
if (!j)c[i][j] = 1;
else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % MOD;
}
int n;
int main() {
//優化輸入
ios::sync_with_stdio(false);
//預處理出來，以后省的再次算了，費一把勁
init();
cin >> n;
//n組詢問
while (n--) {
int a, b;
cin >> a >> b;
printf("%d\n", c[a][b]);
}
return 0;
}
二、預處理法二
例題：https://www.acwing.com/problem/content/888/
上一個辦法是把\(C_a^b\)的值預處理出來了。用的是\(c[N][N]\),\(N\)最大\(2010\).
本題的\(a\),\(b\)都是上限\(10^5\),如果按上題來，就是\(c[10^5][10^5]\),
直接報 \(Memory\ Limit\ Exceeded\). 所以不能直接遞推求出所有解。
適合場景：
1、\(a<=1e5,b<=1e5\)
2、詢問次數:\(\large n<=1e5\)
理論依據：\(\huge C_a^b=\frac{a!}{(a-b)! * b!}\)
舉個栗子：
\(C_3^2=\frac{3 \times 2}{2 \times 1}=3\)
(1) 本題由於怕數據太大，要求結果 \(mod\ (1e9+7)\),這個 \(MOD=1e9+7\)是質數，質數可以直接使用費馬小定理+快速冪求逆元。
(2) 現在要求計算出 \(a!\), \({b!}^{-1}\),還有 \({(a-b)!}^{-1}\)，使用兩個數組來裝(遞推)：
\(fact[i] = i!\ \% \ (1e9+7)\)
\(infact[i] = (i!)^{-1}\ \% \ (1e9+7)\)
所以： \(C_a^b=\frac{a!}{(a-b)! * b!} = (fact[a] * infact[a-b] * infact[b] )\ \% \ (1e9+7)\)
有兩個問題需要突破:
1、需要計算\(a-b\)和\(b\) 在模 \(1e9+7\) 這個乘法世界的逆元
由於\(1e9+7\)是質數，可以利用費馬小定理+快速冪來計算逆元。
2、需要計算\(fact[a]\),即求階乘
求階乘，我們可以在初始化時，從小到大，一路走來一路乘，就可以得到\(a\)的\(fact[a]\)的值，即階乘。只要初始時\(fact[0]=1\)即可,就是零的階乘是\(1\).
C++ 代碼
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 100010; //數據上限
const int MOD = 1e9 + 7; //模值
int fact[N]; //用來保存階乘的值
int infact[N]; //用來保存階乘逆元的值
//快速冪模板
int qmi(int a, int k, int p) {
int res = 1;
while (k) {
if (k & 1) res = (LL) res * a % p;
a = (LL) a * a % p;
k >>= 1;
}
return res;
}
int main() {
//優化輸入
ios::sync_with_stdio(false);
fact[0] = 1; // 0的階乘是1,這是人為的規定。
infact[0] = 1; // 1/1也是1，infact[0]也是1
//對於每一個數字n進行計算
for (int i = 1; i < N; i++) {
// 根據定義求階乘,注意一路要進行MOD運算，防止爆掉
fact[i] = (LL) fact[i - 1] * i % MOD; //強制轉為LL是為了防止越界
// 費馬小定理求i逆元
infact[i] = (LL) infact[i - 1] * qmi(i, MOD - 2, MOD) % MOD;
}
int n;
cin >> n;
while (n--) {
int a, b;
cin >> a >> b;
//公式C(a,b)=a!/(b!*(a-b)!)
printf("%d\n", (LL) fact[a] * infact[b] % MOD * infact[a - b] % MOD);
}
return 0;
}
三、\(Lucas\)公式
例題：https://www.acwing.com/problem/content/889/
1.\(Lucas\)公式:
\(\Huge C_a^b \equiv C_{a\%p}^{b\%p} * C_{\frac{a}{p}}^{\frac{b}{p}} (mod \ p)\)
適用場景：
數據范圍：\(b<=a<=1e18\),\(p<=1e5\),\(p \in prime\),\(n<=20\)
預處理出 \(1...p\) 的階乘和階乘的逆元，用盧卡斯定理進行回答。
2.如何求解\(C_a^b\)
\(C_a^b=\frac{a!}{(a-b)! \times b!}=\frac{a \times (a-1) \times (a-2) \times ...\times(a-b+1) \times (a-b) \times ... \times 1}{(a-b) \times (a-b-1) \times ... \times 1 \times b!}= \frac{a \times (a-1) \times (a-2) \times ...(a-b+1)}{b!}=\frac{a \times (a-1) \times (a-2) \times ...(a-b+1)}{b \times (b-1) \times (b-2) \times ... \times 2 \times 1}\)
根據這個結論，對於\(a\)來講，需要變量\(j\)從\(a\)一直遍歷到\(a-b+1\)。對於\(b\)來講，需要變量\(i\)從\(1\)一起遍歷到\(b\)。而且\(a\)到\(a-b+1\)其實就是\(b\)次，比如\(10\)遍歷到\(8\)，就是\(10,9,8\)三次，即\(b=3\),也就是\(10-3+1=8\)。
而\(i\)也是遍歷\(b\)次，真有太有意思了，它們兩個可以在一個\(b\)次的循環中一起變動！
int C(int a, int b, int p) {
int res = 1;
for (int i = 1, j = a; i <= b; i++, j--) {
//從a到a-b+1是b次，而1-b也是b次，放到一個循環里就行。只不過一個在長大，另一個在減小，這就是上面推導的公式。
res = (LL) res * j % p; // a* (a-1) *(a-2)*....*(a-b+1)
res = (LL) res * qmi(i, p - 2, p) % p; //每次乘i的逆元，而且p是質數，所以可以用費馬小定理快速求逆元
}
return res;
}
3.C++ 代碼
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
/**
* 功能：快速冪模板
* @param a
* @param k
* @param p
* @return
*/
int qmi(int a, int k, int p) {
int res = 1;
while (k) {
if (k & 1) res = (LL) res * a % p;
a = (LL) a * a % p;
k >>= 1;
}
return res;
}
/**
* 功能：組合數模板
* @param a 在a個數中
* @param b 取b個數
* @param p 一個質數，用來取模
* @return 多少種辦法
*/
int C(int a, int b, int p) {
int res = 1;
for (int i = 1, j = a; i <= b; i++, j--) {
//從a到a-b+1是b次，而1-b也是b次，放到一個循環里就行。只不過一個在長大，另一個在減小，這就是上面推導的公式。
res = (LL) res * j % p; // a* (a-1) *(a-2)*....*(a-b+1)
res = (LL) res * qmi(i, p - 2, p) % p; //每次乘i的逆元，而且p是質數，所以可以用費馬小定理快速求逆元
}
return res;
}
/**
* 功能：Lucas公式模板
* @param a
* @param b
* @param p
* @return
*/
int lucas(LL a, LL b, int p) {
if (a < p && b < p) return C(a, b, p);
return (LL) C(a % p, b % p, p) * lucas(a / p, b / p, p) % p; //套用公式，還有個遞歸
}
int n, p;
int main() {
cin >> n;
//n組詢問
while (n--) {
LL a, b;
cin >> a >> b >> p;
//利用lucas公式計算組合數
cout << lucas(a, b, p) << endl;
}
return 0;
}
四、不讓取模怎么辦
https://www.acwing.com/problem/content/890/
1、與前三道題的區別
不再是數據上限的問題了，而是一直不\(mod\),有多大都保留，這無疑可能會爆\(long\ long\),需要使用高精度。
2、公式一【從定義出發的組合數公式】：
(1) \(\large C_a^b=\frac{a \times (a-1) \times ... \times(a-b+1)}{b\times (b-1) \times ... \times 1}=\frac{a \times (a-1) \times ... \times (a-b+1) \times (a-b)!}{ b \times (b-1) \times ... \times 1 \times (a-b)!} =\frac{a!}{b! \times (a-b)!}\)
3、是不是我們利用高精度+上面的組合數公式直接算就行了？
不是的，因為\(yxc\)(大雪菜老師)說，這樣的效率太低，不可以，需要再想一個好辦法。
4、公式二 【算術基本定理】
算術基本定理:\(\large C_a^b=p_1^{\alpha1} \times p_2^{\alpha2} \times p_3^{\alpha3} ... \times p_k^{\alpha k}\)
其中\(p_1\),\(p_2\)...是質數，\(\alpha 1\),\(\alpha 2\),...是指質數\(p_1\),\(p_2\),...的個數。 如果能分解質因數成功的話，那么就可以通過高精度乘法解決掉這個問題。
(1) 那么\(p_1\)和\(p_2\)這些東東怎么求？
思路：因為\(a,b\)都是在\([1,5000]\)之間的數，所以可以通過篩質數的方法，提前獲取到一個質數數組，然后逐個去看，是不是含有這個質數。就能知道有哪些\(p_1,p_2,...\)了。
(2) \(\alpha 1\),\(\alpha 2\)又該怎么求呢？
思路： 求質數\(p\)在\(a!\)中出現的次數辦法：
\(\huge cnt=\lfloor \frac{a}{p} \rfloor + \lfloor \frac{a}{p^2} \rfloor + \lfloor \frac{a}{p^3} \rfloor + ...+ \lfloor \frac{a}{p^k} \rfloor\)
計算\(n\)的階乘中質數\(p\)的個數：
int get(int n, int p) {
int res = 0;
while (n) {
res += n / p;
n /= p;
}
return res;
}
舉個栗子：
$ 5!=1 \times 2 \times 3 \times 4 \times 5$,求\(2\)的個數，
則\(1 \sim 5\)之間含\(2\)的個數就是 \(\lfloor \frac{5}{2} \rfloor\),就是\(2\)個。
則\(1 \sim 5\)之間含\(2^2\)的個數就是 \(\lfloor \frac{5}{2^2} \rfloor\),就是\(1\)個。
\(2^3\)就大於\(5\)，就不再繼續。那么一共的個數就是\(2+1=3\)個。表示在 \(5!\)這個數字中，存在\(3\)個\(2\)，就是\(2\)中有\(1\)個\(2\)，\(4\)中有兩個\(2\)。
5、算法流程
(1)、篩素數，把\(1-5000\)之間的素數篩出來。
(2)、計算\(C_a^b\)中每個已求得素數的個數。
(3)、利用高精度，計算\(C_a^b\)\(=p_1^{\alpha1} \times p_2^{\alpha2} \times p_3^{\alpha3} ... \times p_k^{\alpha k}\)的值。
6、C++代碼
#include <bits/stdc++.h>
using namespace std;
const int N = 5010;
int sum[N]; //每一個質數的次數
int primes[N], cnt;
bool st[N];
//歐拉篩
void get_primes(int n) {
for (int i = 2; i <= n; i++) {
if (!st[i]) primes[cnt++] = i;
for (int j = 0; primes[j] <= n / i; j++) {
st[primes[j] * i] = true;
if (i % primes[j] == 0) break;
}
}
}
//高精度乘法
vector<int> mul(vector<int> a, int b) {
vector<int> c;
int t = 0;
for (int i = 0; i < a.size(); i++) {
t += a[i] * b;
c.push_back(t % 10);
t /= 10;
}
while (t) {
c.push_back(t % 10);
t /= 10;
}
return c;
}
/**
* 功能：n的階乘中包含的質因子p的個數
* @param n n的階乘
* @param p 質因子p
* @return 有多少個
* 參考資料：https://blog.csdn.net/spidy_harker/article/details/88414504
*/
int get(int n, int p) {
int ans = 0;
while (n) { //p^1的個數，p^2的個數,p^3的個數...
ans += n / p;
n /= p;
}
return ans;
}
int main() {
//優化輸入
ios::sync_with_stdio(false);
int a, b;
cin >> a >> b;
//篩2-5000之間質數
get_primes(5000);
//遍歷每一個質數,求每個質數的次數
for (int i = 0; i < cnt; i++) {
int p = primes[i]; //當前質數
//a!中有多少個質數因子p
//減去(a-b)!的多少個質數因子p,
//再減去b!的質數因子p的個數，就是總個數
//sum[i]記錄了每一個質數因子出現的次數
sum[i] = get(a, p) - get(a - b, p) - get(b, p);
}
//利用高精度，把質數因子乘到一起，就是結果了
vector<int> res;
res.push_back(1);
for (int i = 0; i < cnt; i++) //遍歷每一個質數因子
for (int j = 0; j < sum[i]; j++) //開始乘以它的個數
res = mul(res, primes[i]);
//輸出答案
for (int i = res.size() - 1; i >= 0; i--) printf("%d", res[i]);
return 0;
}
7、這段代碼如何理解？
sum[i] = get(a, p) - get(a - b, p) - get(b, p);
看公式一：
\(C_a^b=\frac{a \times (a-1) \times ... \times(a-b+1)}{b\times (b-1) \times ... \times 1}\) \(=\frac{a \times (a-1) \times ... \times (a-b+1) \times (a-b)!}{ b \times (b-1) \times ... \times 1 \times (a-b)!}\) \(=\frac{a!}{b! \times (a-b)!}\)
\(a!\)中質數\(p\)的個數，減去\(b!\)中質數\(p\)的個數，再減去\((a-b)!\)中質數\(p\)的個數，就是公因子消掉的意思。