组合数（C++）

组合数

定义

\qquad 从 $n$ 个不同元素中取出 $m(m\le n)$ 个元素的所有组合的个数，叫做从 $n$ 个不同元素中取出 $m$ 个元素的组合数，记作 $C_n^m$ 或者 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)$。

公式

$$C_n^m=\frac{n(n-1)(n-2)\cdots (n-m+1)}{1\times 2\times 3\times \cdots \times m}=\frac{n!}{m!(n-m)!}$$

\qquad 特别地，$C_n^0=1$.

性质

1. 互补：$C_n^m=C_n^{n-m}$.
2. 组合恒等式：$C_n^m=C_{n-1}^m+C_{n-1}^{m-1}$.

应用

\qquad 针对不同的数据范围应采用不同的算法。

一、885. 求组合数 I - AcWing题库

题目

\qquad 给定 $n$ 组询问，每组询问给定两个整数 $a，b$，请你输出 $C_a^{b}mod(10^9+7)$ 的值。

输入格式

\qquad 第一行包含整数 $n$。

\qquad 接下来 $n$ 行，每行包含一组 $a$ 和 $b$。

输出格式

\qquad 共 $n$ 行，每行输出一个询问的解。

数据范围

\qquad $1\le n\le 1\times 10^4,1\le b\le a\le 2\times 10^3$

分析

\qquad 此题数据量较大，如果对每组数据分别进行计算的话，要算 $2\times 10^7$ 组数据，很可能会超时。但 $a$ 和 $b$ 都较小，如果先预处理出所有可能的 $C_a^b$ 的话，只需要计算 $(2\times 10^3{)}^2/2=2\times 10^6$ 组数据，也就是说 $n$ 组询问里有很多是重复询问的，无需每次都计算一遍。

时间复杂度

\qquad $O(n^2)$

代码

#include <iostream>
using namespace std;

const int N = 2e3 + 5;
const int mod = 1e9 + 7;
int c[N][N];

void init() {				//预处理出所有可能询问的数据
	for (int i = 0; i < N; i++)
		for (int j = 0; j <= i; j++) {
			if (j == 0 || j == i)
				c[i][j] = 1;
			else			//利用组合数恒等式递推计算
				c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;
		}
}

int main() {
	init();
	int n;
	scanf("%d", &n);
	while (n--) {
		int a, b;
		scanf("%d%d", &a, &b);
		printf("%d\n", c[a][b]);
	}
	return 0;
}

二、886. 求组合数 II - AcWing题库

题目

\qquad 给定 $n$ 组询问，每组询问给定两个整数 $a，b$，请你输出 $C_a^{b}mod(10^9+7)$ 的值。

输入格式

\qquad 第一行包含整数 $n$。

\qquad 接下来 $n$ 行，每行包含一组 $a$ 和 $b$。

输出格式

\qquad 共 $n$ 行，每行输出一个询问的解。

数据范围

\qquad $1\le n\le 1\times 10^4,1\le b\le a\le 1\times 10^5$

分析

\qquad 此题与上一题唯一的区别就是 $a$ 和 $b$ 可能的值都变大了，那么我们就不能再使用上一题的思路来求解了。

\qquad 由于 $\frac{a}{b}(\mathrm{mod}p)\ne \frac{a(\mathrm{mod}p)}{b(\mathrm{mod}p)}$，因此我们在用公式求这个 $C_a^b(\mathrm{mod}p)$ 的值的时候就势必要求 $i!$ 的逆元。

\qquad 因此这题我们用公式求解，先预处理出可能用到的 $i!$ 以及 $i!$ 的逆元，分别用 fac[i] 和 infac[i] 数组存储。

\qquad 最终 $C_a^b(\mathrm{mod}p)$ 就可以表示为：
$$\left(fac[a]\ast infac[a-b]\ast infac[b]\right)modp$$
\qquad 还有一点很关键：

\qquad 由于此题的 $mod=10^9+7$，其恰好是一个质数，因此我们利用费马小定理直接用快速幂求逆元即可，即求 infac[i]=fastpow(fac[i],mod-2,mod}.

\qquad 另：这里还可以换种求法，由于 $(n!{)}^{-1}=((n-1)!\times n{)}^{-1}=(n-1){!}^{-1}\times n^{-1}$，所以也可以用infac[i]=infac[i-1]*fastpow(i,mod-2,mod) 这样的方法求解。

时间复杂度

\qquad 预处理 $n$ 个数，每个数都要进行快速幂：$logmod$.

\qquad $O(nlogn)$

代码

#include <iostream>
using namespace std;
typedef long long LL;

const int N = 1e5 + 5;
const int mod = 1e9 + 7;
int fac[N], infac[N];

int fastpow(int base) {						//快速幂
	int res = 1;
	int power = mod - 2;
	while (power) {
		if (power & 1)
			res = (LL)res * base % mod;		//LL不能丢
		base = (LL)base * base % mod;
		power >>= 1;
	}
	return res;
}

void init() {								//预处理出所有可能询问的数据
	fac[0] = infac[0] = 1;
	for (int i = 1; i < N; i++) {
		fac[i] = (LL)fac[i - 1] * i % mod;	//LL不能丢
		infac[i] = fastpow(fac[i]);
	}
}

int main() {
	init();
	int n;
	scanf("%d", &n);
	while (n--) {
		int a, b;
		scanf("%d%d", &a, &b);		//LL不能丢，要及时取模（mod了两次）
		printf("%lld\n", (LL)fac[a] * infac[a - b] % mod * infac[b] % mod);
	}
	return 0;
}

三、887. 求组合数 III - AcWing题库

题目

\qquad 给定 $n$ 组询问，每组询问给定两个整数 $a,b,p$，其中 $p$ 是质数。请你输出 $C_a^b(\mathrm{mod}p)$ 的值。

输入格式

\qquad 第一行包含整数 $n$。

\qquad 接下来 $n$ 行，每行包含一组 $a,b,p$。

输出格式

\qquad 共 $n$ 行，每行输出一个询问的解。

数据范围

\qquad $1\le n\le 20,1\le b\le a\le 1\times 10^{18},1\le p\le 10^5$

分析

\qquad 此题的询问量很小，但是每次都可能询问很大的数字。考虑使用卢卡斯（lucas）定理：

令 $n=sp+q,m=tp+r.(0\le q,r\le p-1)$，

则
$$C_n^m=C_{sp+q}^{tp+r}\equiv C_s^t\cdot C_q^r(\mathrm{mod}p)$$
即
$$C_n^m\equiv C_{n\%p}^{m\%p}\cdot C_{n/p}^{m/p}(\mathrm{mod}p)$$
\qquad 这样处理之后，我们要求的组合数就在 $10^5$ 以内了，可以选择直接用定义求解。但是由于此题 p 不唯一，因此不能沿用第二题的方法预处理出能用到的阶乘和逆元。

代码

#include <iostream>
using namespace std;
typedef long long LL;

int p;

int fastpow(int base) {		//快速幂
	int res = 1;
	int power = p - 2;
	while (power) {
		if (power & 1)
			res = (LL)res * base % p;
		base = (LL)base * base % p;
		power >>= 1;
	}
	return res;
}

int C(int a, int b) {		//用定义求解组合数
	int res = 1;
	for (int i = 1, j = a; i <= b; i++, j--) {
		res = (LL)res * j % p;
		res = (LL)res * fastpow(i) % p;
	}
	return res;
}

int lucas(LL a, LL b) {     //卢卡斯定理的使用
	if (a < p && b < p)	return C(a, b);
	return (LL)C(a % p, b % p) * lucas(a / p, b / p) % p;
}

int main() {
	int n;
	scanf("%d", &n);
	while (n--) {
		LL a, b;
		scanf("%lld%lld%d", &a, &b, &p);
		printf("%d\n", lucas(a, b));
	}
	return 0;
}

四、888. 求组合数 IV - AcWing题库

题目

\qquad 输入 $a,b$，求 $C_a^b$ 的值。

\qquad 注意结果可能很大，需要使用高精度计算。

输入格式

\qquad 共一行，包含两个整数 $a$ 和 $b$。

输出格式

\qquad 共一行，输出 $C_a^b$ 的值。

数据范围

\qquad $1\le b\le a\le 5000$

分析

\qquad 此题与前面截然不同，此题不再需要对答案取模，而是求出精确值。如果按照公式计算，那么需要写一个高精度乘法和一个高精度除法，这是很复杂的。

\qquad 在这种情况下，由于组合数的值一定是一个整数，我们可以对分子分母分别进行质因数分解，然后进行约分，这样就只需要进行高精度乘法运算，而无需再写复杂的高精度除法运算了。

对 n! 进行质因数分解：

\qquad 以 $10!$ 为例，我们想知道将这个数质因数分解后 $2$ 的指数为几（res）。首先我们知道 $10!=1\times 2\times 3\times \cdots \times 10$，其中 $2，4，6，8，10$ 均为 $2$ 的倍数，那么先把这 $5$ 个质因子加上(res += 10/2)；然后我们发现 $4,8$ 中分别有 $2,3$ 个质因数 $2$，而我们只算了一次，通过 res += 10 / 4 就可以把 $4,8$ 分别再加一次；再通过 res += 10 / 8 可以把 $8$ 再加一次。至此，我们用三次运算统计出了 $10!$ 中 $2$ 出现了 $8$ 次。（上述所有除法均为向下取整）

\qquad 总结一下，对于 $n!$ 来说，它分解质因数后质因子 $p$ 的指数为：
$$\lfloor \frac{n}{p^1}\rfloor +\lfloor \frac{n}{p^2}\rfloor +\lfloor \frac{n}{p^3}\rfloor +\cdots +\lfloor \frac{n}{p^k}\rfloor ,(p^k\le n)$$

代码

#include <iostream>
#include <cstring>
#define MOD 10000
using namespace std;

const int N = 1e5 + 5;

int p[N], t;
bool s[N];
int sum[N];				//sum[i]表示第i个质数p[i]出现的次数（指数）

struct HP {				//四位压缩高精度
	int l[2000];
	int len;
	HP() {
		memset(l, 0, sizeof(l));
		len = 0;
	}
	void print() {
        printf("%d",l[len]);
		for (int i = len - 1; i > 0; i--)
			printf("%04d", l[i]);
		printf("\n");
	}
};
HP operator*(const HP &a, const int &b) {		//高精乘低精
	HP c;
	c.len = a.len;
	int x = 0;
	for (int i = 1; i <= c.len; i++) {
		c.l[i] = a.l[i] * b + x;
		x = c.l[i] / MOD;
		c.l[i] %= MOD;
	}
	while (x) {
		c.len++;
		c.l[c.len] = x % MOD;
		x /= MOD;
	}
	return c;
}
void prime(int n) {				//筛n以内的素数
	int i, j;
	for (i = 2; i <= n; i++) {
		if (!s[i])	p[++t] = i;
		for (j = 1; p[j] <= n / i; j++) {
			s[i * p[j]] = 1;
			if (i % p[j] == 0)	break;
		}
	}
}

int get(int n, int p) {			//n!中p这个质数出现了几次
	int res = 0;
	while (n) {
		res += n / p;
		n /= p;
	}
	return res;
}

int main() {
	int a, b;
	cin >> a >> b;
	prime(a);
	for (int i = 1; i <= t; i++) {
		int pr = p[i];
		sum[i] = get(a, pr) - get(b, pr) - get(a - b, pr);		//相当于约分操作
	}
	
	HP res;
    res.l[1] = 1;
    res.len = 1;
	for (int i = 1; i <= t; i++)
		for (int j = 1; j <= sum[i]; j++)
			res = res * p[i];					//做乘法
	res.print();
	return 0;
}