数据结构与算法学习笔记----求组合数

@@ author: 明月清了个风
@@ first publish time: 2025.1.27

ps⭐️一组求组合数的模版题，因为数据范围的不同要用不同的方法进行求解，涉及了很多之前的东西快速幂，逆元，质数，高精度等等，还有新的定理——lucas定理，卡特兰数

Acwing 885. 求组合数 I

[原题链接](885. 求组合数 I - AcWing题库)

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

输入格式

第一行包含整数$n$，

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

输出格式

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

数据范围

$1\le n\le 10000$，

$1\le b\le a\le 2000$

思路

首先我需要知道组合数的基本计算方法，如有一组合数$C_a^b$，那么其计算如下：
$$C_a^b=\frac{a!}{b!(a-b)!}\text{\hspace{1em}(1)}$$
这里引入组合数的一个公式如下：
$$C_a^b=C_{a-1}^b+C_{a-1}^{b-1}\text{\hspace{1em}(2)}$$
我们可以这样理解这个公式：假设我们是要在$a$个苹果中选取$b$个苹果，即$C_a^b$。那么对于某一个苹果而言，只可能有两种情况——被选取和不被选取，这两种情况包含了我们要求的目标的所有方案。当这个苹果被选取了，那么这样的方案数就是$C_{a-1}^{b-1}$，即在剩下的$a-1$个苹果中再选取$b-1$个；当这个苹果没有被选取，其对应方案数就是$C_{a-1}^b$，也就是在剩下的$a-1$个苹果中仍需要选取$b$个。

这样就意味着我们的目标是可以通过之前的组合数计算出来的，也就是可以进行递推，而不用每次都重新计算，时间复杂度也就降低了。

代码

#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>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 ++)if(!j) c[i][j] = 1;else c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mod;
}int main()
{init();int n;cin >> n;while(n --){int a, b;cin >> a >> b;cout << c[a][b] << endl;}return 0;
}

Acwing 886. 求组合数 II

[原题链接](886. 求组合数 II - AcWing题库)

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

输入格式

第一行包含整数$n$，

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

输出格式

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

数据范围

$1\le n\le 10000$，

$1\le b\le a\le 10^5$

思路

注意这题的数据范围$a,b$都到了$10^5$，如果还用上面的预处理两重循环会直接超时，这里我们直接从组合数的计算公式出发看我们要计算的是什么
$$C_a^b\mathrm{mod}(10^9+7)=\frac{a!}{b!(a-b)!}\mathrm{mod}(10^9+7)\text{\hspace{1em}(1)}$$
这样其实就是我们熟悉的内容了，对于阶乘而言我们可以进行预处理，唯一需要注意的就是除法模运算的处理，我们可以使用之前学过的逆元将其转换为一个乘法模运算，需要做的是预处理出分母阶乘的逆元。

逆元用之前的快速幂模版即可求得，具体的看一下快速幂那一篇。

代码

#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>using namespace std;typedef long long LL;const int N = 100010, mod = 1e9 + 7;int n;
int fact[N], infact[N];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()
{cin >> n;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;}while(n --){int a, b;cin >> a >> b;cout << (LL)fact[a]  * infact[a - b] % mod * infact[b] % mod << endl;}return 0;}

Acwing 887. 求组合数 III

[原题链接](887. 求组合数 III - AcWing题库)

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

输入格式

第一行包含整数$n$，

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

输出格式

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

数据范围

$1\le n\le 20$，

$1\le b\le a\le 10^{18}$，

$1\le p\le 10^5$,

思路

到这一题我们可以发现$a$和$b$的数据范围变得非常大，并且模数$p$也不再是一个固定值，因此我们引入一个新的定理——Lucas定理：
$$C_a^b\equiv C_{a\mathrm{mod}p}^{b\mathrm{mod}p}\cdot C_{a/p}^{b/p}(\mathrm{mod}p)$$
也就是我们可以将组合数分解为两个更小的组合数的乘积，并且可以递归分解，直到我们可以处理

这里的证明比较复杂，我们将y总的证明主要分为两点进行说明：

1. 🏫 第一点是对于任意质数$p$，有$(1 + x)^p \equiv 1 + x^p \pmod p $

   这个等式中$p$一定要是质数，这是一个强条件，我们可以通过二项式定理进行证明，我们将$(1+x{)}^p$用二项式定理展开可得：
   $$(1+x{)}^p=\sum _{k=0}^p{C}_p^k\cdot 1^{p-k}\cdot x^k=\sum _{k=0}^p{C}_p^k\cdot x^k=C_p^0+C_p^1\cdot x+C_p^2\cdot x^2+\cdots +C_p^p\cdot x^p$$
   其中，$ C_p^k $ 表示从 $ p $ 个不同元素中取出 $ k $ 个元素的组合数。

   而对于任意质数$p$，当 $ 1 \leq k \leq p-1 $ 时，有 $ p $ 能整除 $ C_p^k $。这是因为：
   $$C_p^k=\frac{p!}{k!(p-k)!}=\frac{p\cdot (p-1)\cdot (p-2)\cdots (p-k+1)}{k!}$$
   由于 $p$ 是质数，且 $k$ 和 $p-k$ 均小于 $p$，因此 $p$ 在分子中作为因子出现，而在分母中没有对应的因子可以消去，所以 $p$ 能整除 $C_p^k$。

   那么由于 $p$ 能整除 $C_p^k$（当 $1\le k\le p-1$ 时），因此在模 $p$ 的意义下，这些项都等于 0：
   $$C_p^k\cdot x^k\equiv 0(\mathrm{mod}p)(1\le k\le p-1)$$
   这意味着在 $(1+x{)}^p$ 的展开式中，除了 $C_p^0$ 和 $C_p^p$ 这两项外，其他所有项在模 $p$ 的意义下都为 0。

   那么我们就可以将$(1+x{)}^p$在模$p$的意义下进行简化：
   $$(1+x{)}^p\equiv C_p^0+C_p^p\cdot x^p\mathrm{mod}p$$
   由于 $C_p^0=1$ 且 $C_p^p=1$（因为从 $ p $ 个元素中取出 $ p $ 个或 0 个元素的组合数都是 1），所以得证：

   $$(1+x{)}^p\equiv 1+x^p\mathrm{mod}p$$

2. 🏫第二点就是将$a$和$b$都写为$p$进制数
   $$a=a_k{p}^k+a_{k-1}{p}^{k-1}+\cdots +a_{1}p+a_0$$

   $$b=b_k{p}^k+b_{k-1}{p}^{k-1}+\cdots +b_{1}p+b_0$$

   则$(1+x{)}^a$可以表示为$(1+x{)}^{a_k{p}^k+a_{k-1}{p}^{k-1}+\cdots +a_{1}p+a_0}$，将指数相加变为底数相乘：
   $$(1+x{)}^{a_k{p}^k}\cdot (1+x{)}^{a_{k-1}{p}^{k-1}}\cdot \dots \cdot (1+x{)}^{a_{1}p}\cdot (1+x{)}^{a_0}$$
   这里可以将$p^k$都先计算，即写成两层指数的形式
   $$(1+x{)}^a=((1+x{)}^{p^k}{)}^{a_k}\cdot ((1+x{)}^{p^{k-1}}{)}^{a_{k-1}}\cdot \cdots \cdot ((1+x{)}^{p^1}{)}^{a_1}\cdot ((1+x{)}^{p^0}{)}^{a_0}$$
   根据上面的第一点我们可以得到
   $$(1+x{)}^a\equiv (1+x^{p^k}{)}^{a_k}\cdot (1+x^{p^{k-1}}{)}^{a_{k-1}}\cdot \cdots \cdot (1+x^{p^1}{)}^{a_1}\cdot (1+x^{p^0}{)}^{a_0}(\mathrm{mod}p)$$
   到这里我们通过对比两遍$x^b$的系数就可以知道$C_a^b\equiv C_{a_k}^{b^k}{C}_{a_{k-1}}^{b^{k-1}}\cdot \cdots \cdot C_{a_1}^{b^1}{C}_{a_0}^{b^0}(\mathrm{mod}p)$

这里给出一个例子吧，可能会好理解一点（应该正确吧）

假设想要计算$C(10,5)\mathrm{mod}3$，其中$a=10$，$b=5$，$p=3$。
将$a$和$b$转换为3进制数：
$10_{10}=101_3$（因为$10=1\cdot 3^2+0\cdot 3^1+1\cdot 3^0$）
$5_{10}=12_3$（因为$5=1\cdot 3^1+2\cdot 3^0$）
应用Lucas定理：
$C(10,5)\equiv C(1,1)\cdot C(0,2)\cdot C(1,0)\mathrm{mod}3$
由于$C(0,2)=0$（因为从0个元素中选取2个元素是不可能的），所以整个表达式的结果为0。
因此，$C(10,5)\mathrm{mod}3=0$。

代码

#include <iostream>
#include <cstring>
#include <cstring>
#include <algorithm>using namespace std;typedef long long LL;int p;int qmi(int a, int k)
{int res = 1;while(k){if(k & 1) res = (LL)res * a % p;a = (LL)a * a % p;k >>= 1;}return res;
}int C(int a, int b)
{if (b > a) return 0;int res = 1;for (int i = 1, j = a; i <= b; i ++, j -- ){res = (LL)res * j % p;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;
}

Acwing 888. 求组合数 IV

[原题链接](887. 求组合数 III - AcWing题库)

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

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

输入格式

共一行，包含$a$,$b$的值。

输出格式

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

数据范围

$1\le b\le a\le 5000$，

思路

这一题的区别是没有了模数，这意味着我们的答案会非常大，因此需要使用高精度来进行计算，我们直接通过计算式进行计算即可$C_a^b=\frac{a!}{(a-b)!b!}$，但是这样我们就需要进行高精度乘法以及高精度除法，因此可以进行一些优化。

优化的思路是进行质因数分解，需要统计出答案数的所有质因数的指数，这里涉及到的一个知识点就是一个数阶乘a!的某个质因子p的指数为$\lfloor \frac{a}{p}\rfloor +\lfloor \frac{a}{p^2}\rfloor +\lfloor \frac{a}{p^3}\rfloor +\cdots$。

这应该挺好理解的，$a!$会涉及到所有的$1\sim a$中的数，那么对于包含了$p$作为质因子的这些数而言，我们首先统计最后能够贡献一次的数的个数，即$\lfloor \frac{a}{p}\rfloor$；但其中也会有能够贡献两次或更多次的数，也就是$p^k(k\ge 2)$为其因子，对于这些数我们在统计一次$p^1$时已经统计了一次了，因此对于包含$p^2$为因子的数，只需再统计一次即可，即加上$\lfloor \frac{a}{p^2}\rfloor$，以此类推。

这里的代码会涉及到之前的筛质数以及高精度乘法，可以先去看一下。

代码

#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <vector>using namespace std;const int N = 5010;int primes[N], cnt;
int sum[N];
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[i * primes[j]] = true;if(i % primes[j] == 0) break;}}
}int get(int n, int p)
{int res = 0;while(n){res += (n / p);n /= p;}return res;
}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;
}int main()
{int a, b;cin >> a  >> b;get_primes(a);for(int i = 0; i < cnt; i ++){int p = primes[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 --) cout << res[i];puts("");return 0;
}

Acwing 889. 满足条件的01序列

[原题链接](889. 满足条件的01序列 - AcWing题库)

给定$n$个$0$和$n$个$1$，它们将按照某种顺序排成长度为$2n$的序列，求他们能排列成的所有序列中，能构满足任意前缀序列中$0$的个数都不少于$1$的个数的序列有多少个。

输出的答案对$10^9+7$取模

输入格式

共一行，包含整数$n$。

输出格式

共一行，包含一个整数，表示答案。

数据范围

$1\le n\le 10^5$，

思路

这一题其实是卡特兰数的一种数学模型，这里就先直接给出他的计算公式了（百度学了一下，应该不止这一种公式，这里就写y总讲的了）
$$h(n)=C_{2n}^n-C_{2n}^{n-1}=\frac{C_{2n}^n}{n+1}$$
我们这里也用y总讲的模型进行讲解——可以将其看为一个路径计数的问题

我们将选$0$看做向左走一步，选$1$看做向上走一步，假设我们一共需要走12步，那么我们一半一半选的终点就是$(6,6)$，并且所有的路径方案数为$C_{2n}^n=C_{12}^6$下图中给出了一条符合题意的路径，大家可以看一下，到达终点前的所有前缀路径中向左走的步数都不小于向上走的步数，图中红线为对称线

要找出的是所有符合这种条件的路径，因此可以反向思考，也就是找到所有不符合条件的路径，在下图中给出了一条不符合题意的橙色路径，图中紫色斜线为另一条关键线，可以看到，当橙色路径第一次碰到紫色斜线时，表明其破坏了题目所要求的条件，我们也由此可以看出若一条路径会碰到或者穿过紫色斜线，其为不符合题意的路径，我们要做的就是找出所有这样的路径。

这里的做法比较巧妙，我们将橙色路径第一次碰到紫色斜线的部分后面进行轴对称，如下图绿色路径所示，所有经过紫色斜线的路径经过这样操作之后，其终点一定为点$(5,7)$，因为在轴对称之前的终点固定是$(6,6)$。

那么就可以得出这样的结论：所有从$(0,0)$出发到达点$(5,7)$的路径都可以由不满足题意的到达$(6,6)$的路径经过上述变换得到，其数量为$C_{2n}^{n-1}=C_{12}^5$。

因此最后要求的路径数目就是$C_{12}^6-C_{12}^5$。

最后实现到代码上我们肯定要进一步化简，也就是$h(n)=C_{2n}^n-C_{2n}^{n-1}=\frac{C_{2n}^n}{n+1}$

看题目给出的数据范围是$10^5$，因此需要使用的是上面组合数第二题的计算方法。

代码

#include <iostream>
#include <algorithm>using namespace std;typedef long long LL;const int N = 100010, mod = 1e9 + 7;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()
{int n;cin >> n;int a = n * 2, b = n;int res = 1;for (int i = a; i > a - b; i -- ) res = (LL)res * i % mod;for (int i = 1; i <= b; i ++ ) res = (LL)res * qmi(i, mod - 2, mod) % mod;res = (LL)res * qmi(n + 1, mod - 2, mod) % mod;cout << res << endl;return 0;
}