组合数公式专题
\(C_a^b\)的多种场景下的求法
一、\(AcWing\) \(885\). 求组合数 \(I\)
理论依据：\(\large C_a^b=C_{a-1}^b+C_{a-1}^{b-1}\)
适合场景：
\(\large a<=2000,b<=2000\)
感性证明：
有\(a\)个苹果，现在需要选出\(b\)个苹果。一共有多少种选法呢？
答：走到第一个苹果面前，面临两个选择：
- 选择它
- 放弃它
如果选择它，将要面对\(a-1\)个苹果中选\(b-1\)个苹果的问题
如果放弃它，将要面对\(a-1\)个苹果中选\(b\)个苹果的问题
根据加法原理，两种选择方法加在一起，就是方法总数，也就是上面的递推式。
#include <bits/stdc++.h>
using namespace std;
const int N = 2010;
const int mod = 1e9 + 7;
int c[N][N];
void init() {
for (int i = 0; i < N; i++) c[i][0] = 1;
for (int i = 1; i < N; i++)
for (int j = 1; j <= i; j++)
c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;
}
int n;
int main() {
init();
scanf("%d", &n);
while (n--) {
int a, b;
scanf("%d %d", &a, &b);
printf("%d\n", c[a][b]);
}
return 0;
}
二、\(AcWing\) \(886\). 求组合数 \(II\)
上一个办法是把\(C_a^b\)的值预处理出来了。用的是\(c[N][N]\),\(N\)最大\(2010\).
本题的\(a\),\(b\)都是上限\(10^5\),如果按上题来，就是\(c[10^5][10^5]\),
直接报 \(Memory\ Limit\ Exceeded\). 所以不能直接递推求出所有解。
适合场景：
\(a<=1e5,b<=1e5\)
理论依据：\(\large C_a^b=\frac{a!}{(a-b)! * b!}\)
举个栗子：
\(C_3^2=\frac{3 \times 2}{2 \times 1}=3\)
-
本题由于怕数据太大，要求结果 \(mod\ (1e9+7)\),这个 \(MOD=1e9+7\)是质数(一般取模的都是质数)，可以直接使用 费马小定理+快速幂 求逆元
-
现在要求计算出 \(a!\), \({b!}^{-1}\), \({(a-b)!}^{-1}\)，使用两个数组来递推：
\(fact[i] = i!\ \% \ MOD\)
\(infact[i] = (i!)^{-1}\ \% \ MOD = ((i-1)!)^{-1}\ * (i)^{-1} \% \ MOD\)
注：根据模运算的性质，\(i!\)的逆元可以表示为\((i-1)!\)的逆元乘以\(i\)的逆元。
所以： \(C_a^b=\frac{a!}{(a-b)! * b!} = (fact[a] * infact[a-b] * infact[b] )\ \% \ MOD\)
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;
}
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;
}
\(AcWing\) \(887\). 求组合数 \(III\)
1.\(Lucas\)公式
\(\large C_a^b \equiv C_{a\%p}^{b\%p} * C_{a \div p}^{b \div p} (mod \ p)\)
适用场景：
数据范围：\(b<=a<=1e18\),\(p<=1e5\),\(p \in prime\)
预处理出 \(1...p\) 的阶乘和阶乘的逆元，用卢卡斯定理进行回答。
注意：
这个\(p\)是小于\(1e5\)的，看到要求模\(1e9+7\)绕行啊！！！不是干哪个的啊！！！适合\(a,b\)很大，\(p\)很小的场景啊！！！
\(Lucas\)定理的理解与实现
实现的时候，只需要对这个东西继续递归调用\(Lucas\)定理即可。也就是，对于\(C_{a÷p}^{b÷p}\)继续调用，最后就会出现一个连乘。
于是我们发现，这个东西又除又模的，很像进制拆解。
没错。
也就是说，这个\(Lucas\)定理就相当于把\(n,m\)变成了\(p\)进制数，然后对\(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) {
if (a < b) return 0;
int down = 1, up = 1;
for (int i = a, j = 1; j <= b; i--, j++) {
up = (LL)up * i % p;
down = (LL)down * j % p;
}
return (LL)up * qmi(down, p - 2, p) % p;
}
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) {
if (a < b) return 0;
int down = 1, up = 1;
for (int i = a, j = 1; j <= b; i--, j++) {
up = (LL)up * i % p;
down = (LL)down * j % p;
}
return (LL)up * qmi(down, p - 2, p) % p;
}
/**
* 功能：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;
}
四、\(AcWing\) \(888\). 求组合数 \(IV\)
不让取模怎么办?
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,...,\alpha k\)又该怎么求呢？
5、公式三【求质数\(p\)在\(a!\)中出现的次数】
\(\large 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\)。
6、算法流程
(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}\)的值。
7、C++代码
#include <bits/stdc++.h>
using namespace std;
const int N = 5010;
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;
}
}
}
//高精乘低精
void mul(int a[], int &al, int b) {
int t = 0;
for (int i = 1; i <= al; i++) {
t += a[i] * b;
a[i] = t % 10;
t /= 10;
}
while (t) {
a[++al] = t % 10;
t /= 10;
}
}
/**
* 功能：n的阶乘中包含的质因子p的个数
* @param n n的阶乘
* @param p 质因子p
* @return 有多少个
*/
int get(int n, int p) {
int cnt = 0;
while (n) { // p^1的个数，p^2的个数,p^3的个数...
cnt += n / p;
n /= p;
}
return cnt;
}
// C(a,b)的结果，高精度保存到c数组，同时，返回c数组的长度len
void C(int a, int b, int c[], int &cl) {
//高精度的基底，乘法的基数是1
c[1] = 1;
cl = 1; //由于高精度数组中只有一位，是1，所以长度也是1
for (int i = 0; i < cnt; i++) { //枚举区间内所有质数
int p = primes[i];
/*
C(a,b)=a!/(b! * (a-b)!)
a!中有多少个质数因子p
减去(a-b)!的多少个质数因子p,
再减去b!的质数因子p的个数，就是总个数
s记录了p这个质数因子出现的次数
*/
int s = get(a, p) - get(b, p) - get(a - b, p);
while (s--) mul(c, cl, p); // 不断的乘p,结果保存到数组c中。len将带回c的有效长度
}
}
int a, b;
int c[N], cl;
int main() {
cin >> a >> b;
//筛质数
get_primes(N - 1);
C(a, b, c, cl);
//输出答案
for (int i = cl; i >= 1; i--) printf("%d", c[i]);
return 0;
}
8、这段代码如何理解？
s = 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\)的个数，就是公因子消掉的意思。
