组合计数
求组合数 I —— 预处理组合数
给定 \(n\) 组询问，每组询问给定两个整数 \(a\)，\(b\)，请你输出 \(C_a^b\bmod(10^9+7)\) 的值。
数据范围：
$ 1 \le n \le 10000$
$ 1 \le b \le a \le 2000$
如果直接暴力做的话，每次计算一个组合数，最坏需要2000次运算（因为a与b的范围是2000），而最多有10000次询问，所以时间复杂度达到2千万（ \(2000 * 10000 = 2 * 10^7\) ），会超时，所以不能这么算。
我们注意到 \(a\) 和 \(b\) 的范围是2000，所以不同的 \(a\) 和 \(b\) 的方案数有 \(2000^2\) 组，总共只有四百万对（ \(2000^2 = 4 * 10^6\) ），所以我们可以先预处理出所有方案的情况，利用以下递推式：
\(C_a^b\) 表示从 \(a\) 个苹果中选 \(b\) 个的方案数，假设 \(a\) 个苹果中有一个红色苹果，我们可以将方案分为两大类：
- 第一类的方案中选择的 \(b\) 个苹果包含红色苹果，则这种情况下的方案数就是除了红色苹果之外，从剩下的苹果中再选择 \(b-1\) 个苹果： \(C_{a-1}^{b-1}\)
- 第二类的方案中选择的 \(b\) 个苹果不包含红色苹果，则这种情况下的方案数就是在除去红色苹果的剩余部分选择 \(b\) 个苹果：\(C_{a-1}^b\)
时间复杂度：\(O(N^2)\)
#include <iostream>
using namespace std;
const int N = 2010, mod = 1e9 + 7;
int c[N][N];
// 预处理出所有方案
void init()
{
for (int i = 0; i < N; i ++ )
for (int j = 0; j <= i; j ++ )
// 处理边界情况，从i个苹果里选0个苹果的方案数是1
if (!j) 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;
}
求组合数 II —— 预处理阶乘
给定 \(n\) 组询问，每组询问给定两个整数 \(a\)，\(b\)，请你输出 \(C_a^b\bmod(10^9+7)\) 的值。
数据范围：
$ 1 \le n \le 10000$
$ 1 \le b \le a \le 10^5$
此时如果用上一道题的预处理方法，则由于 \(a\) 和 \(b\) 的范围是十万，不同的 \(a\) 和 \(b\) 的方案数有一百亿组（ \(10^{10}\) ），超时QAQ。所以我们需要换一种预处理方式。我们知道组合数公式为：
\[C_n^m = \frac{n!}{m!(n-m)!} \]所以可以将公式中的阶乘预处理出来，\(fact[i] = i! \bmod (10^9+7)\)
又因为
可以通过求逆元 \(x\)（ \(\frac{a}{b} \bmod p = a x \bmod p\) ），将除法转换为乘法，再预处理出 \(infact[i] = (i!)^{-1} \bmod (10^9+7)\)，则
\[C_n^m = fact[n] \times infact[n-m] \times infact[m] \]本题中模数是质数，符合费马小定理的条件，可以使用快速幂求逆元。
快速幂的时间复杂度为\(O(\log k)\)，\(k\) 是幂次数（我们将幂次分为 \(2^0, 2^1,2^2 \cdots\)）
时间复杂度： \(O(N \log N)\)
#include <iostream>
using namespace std;
typedef long long LL;
const int N = 100010;
const int mod = 1e9 + 7;
int n;
LL fact[N], infact[N];
// 快速幂，返回a^k % p
int qmi(int a, int k, int p)
{
int res = 1;
while (k) {
if (k & 1) res = (LL)res * a % p;
k >>= 1;
a = (LL)a * a % p;
}
return res;
}
int main()
{
// 预处理阶乘数组
fact[0] = infact[0] = 1;
for (int i = 1; i < N; i ++ ) {
fact[i] = (LL)fact[i - 1] * i % mod;
infact[i] = (LL)infact[i - 1] * qmi(i, mod - 2, mod) % mod;
}
scanf("%d", &n);
while (n -- ) {
int a, b;
scanf("%d%d", &a, &b);
// 两个10^9相乘不会溢出long long，三个相乘有可能溢出，所以中间取模一次
printf("%d\n", fact[a] * infact[b] % mod * infact[a - b] % mod);
}
return 0;
}
求组合数 III —— 卢卡斯定理（Lucas）
给定 n 组询问，每组询问给定三个整数 \(a,b,p\)，其中 \(p\) 是质数，请你输出 \(C_b^a\bmod p\) 的值。
数据范围：
$ 1 \le n \le 20$
$ 1 \le b \le a \le 10^{18}$
$ 1 \le p \le 10^5$
Lucas 定理
若 \(p\) 是质数，则对于任意整数 \(1 \le m \le n\) ，有:
也就是把 \(n\) 和 \(m\) 表示成 \(p\) 进制数，对 \(p\) 进制下的每一位分别计算组合数，最后再乘起来。证明一般需要用到生成函数知识，此处仅简略证明。
证明：
把 \(n\) 和 \(m\) 表示成 \(p\) 进制数
由二项式定理 \((a+b)^n = \sum\limits_{k=0}^{n} C_n^ka^kb^{n-k}\) 可以得到：
\[\begin{align*} (1+x)^p & = C_p^0 \cdot x^0 \cdot 1^p + C_p^1 \cdot x^1 \cdot 1^{p-1} + \dots + C_p^p \cdot x^p \cdot 1^0 \\ & = C_p^0 \cdot 1 + C_p^1 \cdot x + \dots + C_p^p x^p \end{align*} \]上式第1项到第n-1项分子都有 \(p\) ，对 \(p\) 取模都得零，于是：
\[(1+x)^p \equiv 1+x^p \pmod p \]将 \(n\) 分解为 \(p\) 进制数，我们得到：
\[\begin{align*} (1+x)^n & = (1+x)^{n_0p^0 + n_1p^1 + n_2p^2 + \dots + n_kp^k} \\ & = (1+x)^{n_0} \cdot ((1+x)^p)^{n_1} + ((1+x)^{p^2})^{n_2} + \dots + ((1+x)^{p^k})^{n_k} \\ & \equiv (1+x)^{n_0} \cdot (1+x^{p})^{n_1} \cdot (1+x^{p^2})^{n_2} + \dots + (1+x^{p_k})^{n_k} \end{align*} \]二项式定理展开，上式左边 \(x^m\) 的系数是 \(C_n^m\) ，右边的系数是 \(C_{n_k}^{m_k} \cdot C_{n_{k-1}}^{m_{k-1}} \dots C_{n_0}^{m_0}\)，则
\[C_n^m \equiv C_{n_k}^{m_k} \cdot C_{n_{k-1}}^{m_{k-1}} \dots C_{n_0}^{m_0} \pmod p \]卢卡斯定理就是上式的递归形式。
证毕
将 \(n\) 和 \(m\) 表示成 \(p\) 进制数的时间复杂度是 \(O(\log_p 10^{18})\) ，利用快速幂从定义出发（阶乘除法）计算一个 \(10^5\) 大小的组合数的时间复杂度是 \(O(p \log p)\)，而当 \(p\) 取极限值时，\(\log_p 10^{18}\) 很小，所以总的时间复杂度约为：\(O(p \log p)\)
#include <iostream>
using namespace std;
typedef long long LL;
int p;
// 快速幂，返回a^k % p，p在本题是全局变量，所以不需要作为参数引入
LL qmi(int a, int k)
{
LL res = 1;
while (k) {
if (k & 1) res = (LL)res * a % p;
k >>= 1;
a = (LL)a * a % p;
}
return res;
}
// 返回组合数C_a^b = a! / b! * (a-b)!
int C(int a, int b)
{
int res = 1;
for (int i = 1, j = a; i <= b; i ++, j -- ) {
// 组合数中a!与(a-b)!约分，则分子只需累乘a到a-b+1
res = (LL)res * j % p;
// 分母累积b! 等价于乘以对应的逆元
res = (LL)res * qmi(i, p - 2) % 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;
cin >> n;
while (n -- ) {
LL a, b;
cin >> a >> b >> p;
cout << lucas(a, b) << endl;
}
return 0;
}
求组合数 IV —— 高精度
输入 \(a,b\)，求 \(C_a^b\) 的值。
注意结果可能很大，需要使用高精度计算。
数据范围：
\(1 \le b \le a \le 5000\)
如果从定义出发，我们需要使用高精度乘法和高精度除法，代码实现比较复杂，所以我们首先对组合数分解质因数，从而只需用高精度乘法即可计算结果。
\[C_a^b = \frac{a!}{b!(a-b)!} = p_1^{\alpha_1} \cdot p_2^{\alpha_2} \cdot p_3^{\alpha_3} + \dots \]而 \(a!\) 内包含 \(p\) 的次数为：
\[\lfloor \frac{a}{p} \rfloor + \lfloor \frac{a}{p^2} \rfloor + \lfloor \frac{a}{p^3} \rfloor + \dots \]\(\lfloor \frac{a}{p} \rfloor\) 表示 \(1\) ～ \(a\) 中 \(p\) 的倍数有多少个
\(\lfloor \frac{a}{p^2} \rfloor\) 表示 \(1\) ～ \(a\) 中 \(p^2\) 的倍数有多少个
以此类推...
举例说明，\(1,2,3,4,5,6,7,8\) 中，在 \(2,4,6,8\) 中均出现了4次 \(2\) ，在 \(4,8\) 中出现了2次 \(2^2\) ，在 \(8\) 中出现了1次 \(2^3\)，于是 \(8!\) 中共包含7个 \(2\)。
该步骤时间复杂度是 \(O(\log p)\)
算法思路：
1. 筛质数，将1～5000内的质数筛出来
2. 求每个质数的次数
3. 用高精度乘法把所有的质因子乘到一起
#include <iostream>
#include <vector>
using namespace std;
const int N = 5010;
int primes[N], cnt; // 存储筛选出的质数与对应的数组下标
bool st[N]; // 标记是质数还是合数
int sum[N]; // 存储阶乘中包含的每个质数的次数
// 线性筛质数
void get_primes(int n)
{
// 质数从2开始
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;
}
}
}
// 返回a!里面包含的p的次数
int get(int a, int p)
{
int res = 0;
while (a) {
res += a / p; // a/p + a/p^2 + a/p^3 + ...
a /= p; // 当a < p^k时, a = a / p^k = 0, 跳出循环
}
return res;
}
// 高精度乘法模板
vector<int> mul(vector<int> a, int &b)
{
vector<int> c;
int t = 0;
for (int i = 0; i < a.size() || t; i ++ ) {
if (i < a.size()) t += a[i] * b;
c.push_back(t % 10);
t /= 10;
}
while (c.size() > 1 && c.back() == 0) c.pop_back();
return c;
}
int main()
{
int a, b;
cin >> a >> b;
// 筛出1～N范围内的质数
get_primes(a);
// 求每个质数的次数
for (int i = 0; i < cnt; i ++ ) {
int p = primes[i];
sum[i] = get(a, p) - get(b, p) - get(a - 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 -- ) cout << res[i];
puts("");
return 0;
}
满足条件的01序列 —— 卡特兰数
01排列的方案数与从(0,0)走到(6,6)的路径数一致。任意前缀中0的个数都要小于1的个数，路径上的点的横纵坐标需满足不经过下图红色线上的点。
总共有12步，其中6步向上走，所以从(0,0)走到(6,6)的所有路径数是 \(C_{12}^6\) 。
任何一条从(0,0)走到(6,6)并且经过红边的路径数量（从第一个与红边相交的点做轴对称，可以对应一条从(0,0)走到(5,7)的路径），从而等于从(0,0)走到(5,7)的路径数量 \(C_{12}^5\) 。
所以满足条件的路径数为：\(C_{12}^6 - C_{12}^5\)
通用形式为：
Catalan数列
给定 \(n\) 个 \(0\) 和 \(n\) 个 \(1\)，它们将按照某种顺序排成长度为 \(2n\) 的序列，求它们能排列成的所有序列中，能够满足任意前缀序列中 \(0\) 的个数都不少于 \(1\) 的个数的序列的数量为：
#include <iostream>
using namespace std;
typedef long long LL;
const int mod = 1e9 + 7;
// 模数为质数，所以可以根据费马小定理，用快速幂求逆元
int qmi(int a, int k, int p)
{
int res = 1;
while (k) {
if (k & 1) res = (LL)res * a % mod;
a = (LL)a * a % mod;
k >>= 1;
}
return res;
}
int main()
{
int n;
cin >> n;
int a = 2 * n, b = n;
int res = 1;
// (2n)! / n!
for (int i = a; i > a - b; i -- ) res = (LL)res * i % mod;
// 1 / n!
for (int i = 1; i <= b; i ++ ) res = (LL)res * qmi(i, mod - 2, mod) % mod;
// 1 / n + 1
res = (LL)res * qmi(n + 1, mod - 2, mod) % mod;
cout << res << endl;
return 0;
}