如何求组合数(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小~
if (!j)c[i][j] = 1; //特判，从i个苹果里选择0个苹果，只有一种方案，就是，啥也不选。也可以理解为把二维数组的第一列全都配置为1，做为递推的起始值
else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % MOD;
}
int n;
int main() {
//预处理出来，以后省的再次算了，费一把劲
init();
cin >> n;
//n组询问
while (n--) {
int a, b;
cin >> a >> b;
printf("%d
", 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)
(1)公式:(huge C_n^m=frac{n!}{(n-m)! * m!}=frac{n(n-1)(n-2).....(n-m+1)}{m(m-1)(m-2) imes .... imes 1})
举个栗子：
(C_3^2=frac{3 imes 2}{2 imes 1}=3)
(2) 同时，本题由于怕数据太大，要求结果 (mod (1e9+7)),这个 (1e9+7)是质数，质数可以直接使用费马小定理求逆元。
(3) 现在要求计算出 (n!), ({m!}^{-1}),还有 ({(n-m)!}^{-1})，使用两个数组来装：
(fact[i] = i! \% (1e9+7))
(infact[i] = (i!)^{-1} \% (1e9+7))
所以： (C_n^m=frac{n!}{(n-m)! * m!} = (fact[n] * infact[n-m] * infact[m] ) \% (1e9+7))
有两个问题需要突破:
1、需要计算(n-m)和(m) 在模 (1e9+7) 这个乘法世界的逆元
由于(1e9+7)是质数，可以利用费马小定理+快速幂来计算逆元。
2、需要计算(fact[n]),即求阶乘
求阶乘，我们可以在初始化时，从小到大，一路走来一路乘，就可以得到(n)的(fact[n])的值，即阶乘。只要初始时(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() {
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; //这里使用了乘法MOD的分配性质
}
int n;
cin >> n;
while (n--) {
int a, b;
cin >> a >> b;
printf("%d
", (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)! imes b!}=frac{a imes (a-1) imes (a-2) imes ... imes(a-b+1) imes (a-b) imes ... imes 1}{(a-b) imes (a-b-1) imes ... imes 1 imes b!}= frac{a imes (a-1) imes (a-2) imes ...(a-b+1)}{b!}=frac{a imes (a-1) imes (a-2) imes ...(a-b+1)}{b imes (b-1) imes (b-2) imes ... imes 2 imes 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 1e9+7)这个条件，就是一直不(mod),有多大都保留，这无疑可能会爆(long long),所以需要使用高精度。
2、公式一【从定义出发的组合数公式】：
(1) (large C_a^b=frac{a imes (a-1) imes ... imes(a-b+1)}{b imes (b-1) imes ... imes 1}=frac{a imes (a-1) imes ... imes (a-b+1) imes (a-b)!}{ b imes (b-1) imes ... imes 1 imes (a-b)!} =frac{a!}{b! imes (a-b)!})
3、是不是我们利用高精度+上面的组合数公式直接算就行了？
不是的，因为(yxc)(大雪菜老师)说，这样的效率太低，不可以，需要再想一个好办法。
4、公式二 【分解质因数】
分解质因数:(large C_a^b=p_1^{alpha1} imes p_2^{alpha2} imes p_3^{alpha3} ... imes 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])之间的数，所以我们可以通过筛质数的方法，提前获取到一个质数数组，然后逐个去看，是不是含有这个质数。
(2) (alpha 1),(alpha 2)又该怎么求呢？
思路： 求质数p在a中出现的次数办法：
(cnt=lfloor frac{a}{p}
floor + lfloor frac{a}{p^2}
floor + lfloor frac{a}{p^3}
floor + ...lfloor frac{a}{p^k}
floor)
第一个表示 (frac{a}{p})的下取整，表示(a)中有几个(p),就是(p)的倍数。
第二个表示 表示(a)中有几个(p^2),就是(p^2)的倍数。
以此类推..不重不漏，好神奇！
计算(n)的阶乘中质数(p)的个数：
int get(int n, int p) {
int res = 0;
while (n) {
res += n / p;
n /= p;
}
return res;
}
举个栗子：
$ 5!=1 imes 2 imes 3 imes 4 imes 5$,求(2)的个数，
则(1~5)之间含(2)的个数就是 (lfloor frac{5}{2}
floor),就是(2)个。
则(1~5)之间含(2^2)的个数就是 (lfloor frac{5}{2^2}
floor),就是(1)个。
(2^3)就大于(5)，就不再继续。那么一共的个数就是(2+1=3)个。表示在 (5!)这个数字中，存在(3)个(2)，就是(2)中有(1)个(2)，(4)中有两个(2)。
要是还不太理解，就看看这里吧：
https://blog.csdn.net/spidy_harker/article/details/88414504
5、算法流程
(1)、筛素数，把(1-5000)之间的素数筛出来。
(2)、计算(C_a^b)中每个已求得素数的个数。
(3)、利用高精度，计算(C_a^b)(=p_1^{alpha1} imes p_2^{alpha2} imes p_3^{alpha3} ... imes 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的个数
//也可以理解为是一个模板
//https://blog.csdn.net/spidy_harker/article/details/88414504
int get(int n, int p) {
int res = 0;
while (n) {
res += n / p;
n /= p;
}
return res;
}
int main() {
int a, b;
cin >> a >> b;
//筛2-a之间素数
get_primes(a);
//遍历每一个质数,求每个质数的次数
for (int i = 0; i < cnt; i++) {
int p = primes[i]; //当前质数
//a!中有多少个质数因子p
//减去(a-b)!的多少个质数因子p,
//再减去b!的质数因子p的个数，就是总个数
sum[i] = get(a, p) - get(a - b, p) - get(b, p);
//sum[i]记录了每一个质数因子出现的次数
}
//利用高精度，把质数因子乘到一起，就是结果了
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 imes (a-1) imes ... imes(a-b+1)}{b imes (b-1) imes ... imes 1}) (=frac{a imes (a-1) imes ... imes (a-b+1) imes (a-b)!}{ b imes (b-1) imes ... imes 1 imes (a-b)!}) (=frac{a!}{b! imes (a-b)!})
(a!)中质数(p)的个数，减去(b!)中质数(p)的个数，再减去((a-b)!)中质数(p)的个数，就是公因子消掉的意思。