数据结构与算法学习笔记----求组合数
@@ author: 明月清了个风
@@ first publish time: 2025.1.27ps⭐️一组求组合数的模版题，因为数据范围的不同要用不同的方法进行求解，涉及了很多之前的东西快速幂，逆元，质数，高精度等等，还有新的定理——lucas定理，卡特兰数
Acwing 885. 求组合数 I
[原题链接](885. 求组合数 I - AcWing题库)
给定 n n n组询问，每组询问给定两个整数 a a a， b b b，请你输出 C a b m o d ( 1 0 9 + 7 ) C_{a}^{b} \bmod (10^9 + 7) Cabmod(109+7)的值。
输入格式
第一行包含整数 n n n，
接下来 n n n行，每行包含一组 a a a和 b b b。
输出格式
共 n n n行，每行输出一个询问的解。
数据范围
1 ≤ n ≤ 10000 1 \le n \le 10000 1≤n≤10000，
1 ≤ b ≤ a ≤ 2000 1 \le b \le a \le 2000 1≤b≤a≤2000
思路
首先我需要知道组合数的基本计算方法，如有一组合数
C
a
b
C_{a}^{b}
Cab，那么其计算如下：
C
a
b
=
a
!
b
!
(
a
−
b
)
!
(1)
C_{a}^{b} = \frac{a!}{b!(a - b)!} \tag{1}
Cab=b!(a−b)!a!(1)
这里引入组合数的一个公式如下：
C
a
b
=
C
a
−
1
b
+
C
a
−
1
b
−
1
(2)
C_{a}^{b} = C_{a - 1}^{b} + C_{a - 1}^{b - 1} \tag{2}
Cab=Ca−1b+Ca−1b−1(2)
我们可以这样理解这个公式：假设我们是要在
a
a
a个苹果中选取
b
b
b个苹果，即
C
a
b
C_{a}^{b}
Cab。那么对于某一个苹果而言，只可能有两种情况——被选取和不被选取，这两种情况包含了我们要求的目标的所有方案。当这个苹果被选取了，那么这样的方案数就是
C
a
−
1
b
−
1
C_{a - 1}^{b - 1}
Ca−1b−1，即在剩下的
a
−
1
a - 1
a−1个苹果中再选取
b
−
1
b - 1
b−1个；当这个苹果没有被选取，其对应方案数就是
C
a
−
1
b
C_{a - 1}^{b}
Ca−1b，也就是在剩下的
a
−
1
a - 1
a−1个苹果中仍需要选取
b
b
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 n n组询问，每组询问给定两个整数 a a a， b b b，请你输出 C a b m o d ( 1 0 9 + 7 ) C_{a}^{b} \bmod (10^9 + 7) Cabmod(109+7)的值。
输入格式
第一行包含整数 n n n，
接下来 n n n行，每行包含一组 a a a和 b b b。
输出格式
共 n n n行，每行输出一个询问的解。
数据范围
1 ≤ n ≤ 10000 1 \le n \le 10000 1≤n≤10000，
1 ≤ b ≤ a ≤ 1 0 5 1 \le b \le a \le 10^5 1≤b≤a≤105
思路
注意这题的数据范围
a
,
b
a,b
a,b都到了
1
0
5
10^5
105，如果还用上面的预处理两重循环会直接超时，这里我们直接从组合数的计算公式出发看我们要计算的是什么
C
a
b
m
o
d
(
1
0
9
+
7
)
=
a
!
b
!
(
a
−
b
)
!
m
o
d
(
1
0
9
+
7
)
(1)
C_{a}^{b} \bmod (10^9 + 7) = \frac{a!}{b!(a - b)!} \bmod (10^9 + 7) \tag{1}
Cabmod(109+7)=b!(a−b)!a!mod(109+7)(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 n n组询问，每组询问给定两个整数 a a a， b b b， p p p，其中 p p p是质数，请你输出 C a b m o d p C_{a}^{b} \bmod p Cabmodp的值。
输入格式
第一行包含整数 n n n，
接下来 n n n行，每行包含一组 a a a, b b b, p p p。
输出格式
共 n n n行，每行输出一个询问的解。
数据范围
1 ≤ n ≤ 20 1 \le n \le 20 1≤n≤20，
1 ≤ b ≤ a ≤ 1 0 18 1 \le b \le a \le 10^{18} 1≤b≤a≤1018，
1 ≤ p ≤ 1 0 5 1 \le p \le 10^5 1≤p≤105,
思路
到这一题我们可以发现
a
a
a和
b
b
b的数据范围变得非常大，并且模数
p
p
p也不再是一个固定值，因此我们引入一个新的定理——Lucas定理：
C
a
b
≡
C
a
m
o
d
p
b
m
o
d
p
⋅
C
a
/
p
b
/
p
(
m
o
d
p
)
C_{a}^{b} \equiv C_{a \mod p}^{b \mod p} \cdot C_{a / p}^{b / p} \pmod p
Cab≡Camodpbmodp⋅Ca/pb/p(modp)
也就是我们可以将组合数分解为两个更小的组合数的乘积，并且可以递归分解，直到我们可以处理
这里的证明比较复杂，我们将y总的证明主要分为两点进行说明：
-
🏫 第一点是对于任意质数 p p p，有$(1 + x)^p \equiv 1 + x^p \pmod p $
这个等式中 p p p一定要是质数，这是一个强条件，我们可以通过二项式定理进行证明，我们将 ( 1 + x ) p (1 + x)^p (1+x)p用二项式定理展开可得：
( 1 + x ) p = ∑ k = 0 p C p k ⋅ 1 p − k ⋅ x k = ∑ k = 0 p C p k ⋅ x k = C p 0 + C p 1 ⋅ x + C p 2 ⋅ x 2 + ⋯ + C p p ⋅ 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 (1+x)p=k=0∑pCpk⋅1p−k⋅xk=k=0∑pCpk⋅xk=Cp0+Cp1⋅x+Cp2⋅x2+⋯+Cpp⋅xp
其中，$ C_p^k $ 表示从 $ p $ 个不同元素中取出 $ k $ 个元素的组合数。而对于任意质数 p p p，当 $ 1 \leq k \leq p-1 $ 时，有 $ p $ 能整除 $ C_p^k $。这是因为：
C p k = p ! k ! ( p − k ) ! = p ⋅ ( p − 1 ) ⋅ ( p − 2 ) ⋯ ( p − k + 1 ) k ! C_p^k = \frac{p!}{k!(p-k)!} = \frac{p \cdot (p-1) \cdot (p-2) \cdots (p-k+1)}{k!} Cpk=k!(p−k)!p!=k!p⋅(p−1)⋅(p−2)⋯(p−k+1)
由于 p p p 是质数，且 k k k 和 p − k p-k p−k 均小于 p p p，因此 p p p 在分子中作为因子出现，而在分母中没有对应的因子可以消去，所以 p p p 能整除 C p k C_p^k Cpk。那么由于 p p p 能整除 C p k C_p^k Cpk（当 1 ≤ k ≤ p − 1 1 \leq k \leq p-1 1≤k≤p−1 时），因此在模 p p p 的意义下，这些项都等于 0：
C p k ⋅ x k ≡ 0 ( m o d p ) ( 1 ≤ k ≤ p − 1 ) C_p^k \cdot x^k \equiv 0 \pmod p \quad (1 \leq k \leq p-1) Cpk⋅xk≡0(modp)(1≤k≤p−1)
这意味着在 ( 1 + x ) p (1+x)^p (1+x)p 的展开式中，除了 C p 0 C_p^0 Cp0 和 C p p C_p^p Cpp 这两项外，其他所有项在模 p p p 的意义下都为 0。那么我们就可以将 ( 1 + x ) p (1 + x)^p (1+x)p在模 p p p的意义下进行简化：
( 1 + x ) p ≡ C p 0 + C p p ⋅ x p m o d p (1+x)^p \equiv C_p^0 + C_p^p \cdot x^p \mod p (1+x)p≡Cp0+Cpp⋅xpmodp
由于 C p 0 = 1 C_p^0 = 1 Cp0=1 且 C p p = 1 C_p^p = 1 Cpp=1（因为从 $ p $ 个元素中取出 $ p $ 个或 0 个元素的组合数都是 1），所以得证：( 1 + x ) p ≡ 1 + x p m o d p (1+x)^p \equiv 1 + x^p \mod p (1+x)p≡1+xpmodp
-
🏫第二点就是将 a a a和 b b b都写为 p p p进制数
a = a k p k + a k − 1 p k − 1 + ⋯ + a 1 p + a 0 a = a_k p^k + a_{k-1} p^{k-1} + \cdots + a_1 p + a_0 a=akpk+ak−1pk−1+⋯+a1p+a0b = b k p k + b k − 1 p k − 1 + ⋯ + b 1 p + b 0 b = b_k p^k + b_{k-1} p^{k-1} + \cdots + b_1 p + b_0 b=bkpk+bk−1pk−1+⋯+b1p+b0
则 ( 1 + x ) a (1 + x)^a (1+x)a可以表示为 ( 1 + x ) a k p k + a k − 1 p k − 1 + ⋯ + a 1 p + a 0 (1+x)^{a_k p^k + a_{k-1} p^{k-1} + \cdots + a_1 p + a_0} (1+x)akpk+ak−1pk−1+⋯+a1p+a0，将指数相加变为底数相乘：
( 1 + x ) a k p k ⋅ ( 1 + x ) a k − 1 p k − 1 ⋅ … ⋅ ( 1 + x ) a 1 p ⋅ ( 1 + x ) a 0 (1+x)^{a_k p^k} \cdot (1+x)^{a_{k-1} p^{k-1}} \cdot \ldots \cdot (1+x)^{a_1 p} \cdot (1+x)^{a_0} (1+x)akpk⋅(1+x)ak−1pk−1⋅…⋅(1+x)a1p⋅(1+x)a0
这里可以将 p k p^{k} pk都先计算，即写成两层指数的形式
( 1 + x ) a = ( ( 1 + x ) p k ) a k ⋅ ( ( 1 + x ) p k − 1 ) a k − 1 ⋅ ⋯ ⋅ ( ( 1 + x ) p 1 ) a 1 ⋅ ( ( 1 + x ) p 0 ) a 0 (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=((1+x)pk)ak⋅((1+x)pk−1)ak−1⋅⋯⋅((1+x)p1)a1⋅((1+x)p0)a0
根据上面的第一点我们可以得到
( 1 + x ) a ≡ ( 1 + x p k ) a k ⋅ ( 1 + x p k − 1 ) a k − 1 ⋅ ⋯ ⋅ ( 1 + x p 1 ) a 1 ⋅ ( 1 + x p 0 ) a 0 ( m o d p ) (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} \pmod p (1+x)a≡(1+xpk)ak⋅(1+xpk−1)ak−1⋅⋯⋅(1+xp1)a1⋅(1+xp0)a0(modp)
到这里我们通过对比两遍 x b x^{b} xb的系数就可以知道 C a b ≡ C a k b k C a k − 1 b k − 1 ⋅ ⋯ ⋅ C a 1 b 1 C a 0 b 0 ( m o d p ) 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} \pmod p Cab≡CakbkCak−1bk−1⋅⋯⋅Ca1b1Ca0b0(modp)
这里给出一个例子吧，可能会好理解一点（应该正确吧）
假设想要计算 C ( 10 , 5 ) m o d 3 C(10, 5) \mod 3 C(10,5)mod3，其中 a = 10 a = 10 a=10， b = 5 b = 5 b=5， p = 3 p = 3 p=3。
将 a a a和 b b b转换为3进制数：
1 0 10 = 10 1 3 10_{10} = 101_3 1010=1013（因为 10 = 1 ⋅ 3 2 + 0 ⋅ 3 1 + 1 ⋅ 3 0 10 = 1 \cdot 3^2 + 0 \cdot 3^1 + 1 \cdot 3^0 10=1⋅32+0⋅31+1⋅30）
5 10 = 1 2 3 5_{10} = 12_3 510=123（因为 5 = 1 ⋅ 3 1 + 2 ⋅ 3 0 5 = 1 \cdot 3^1 + 2 \cdot 3^0 5=1⋅31+2⋅30）
应用Lucas定理：
C ( 10 , 5 ) ≡ C ( 1 , 1 ) ⋅ C ( 0 , 2 ) ⋅ C ( 1 , 0 ) m o d 3 C(10, 5) \equiv C(1, 1) \cdot C(0, 2) \cdot C(1, 0) \mod 3 C(10,5)≡C(1,1)⋅C(0,2)⋅C(1,0)mod3
由于 C ( 0 , 2 ) = 0 C(0, 2) = 0 C(0,2)=0（因为从0个元素中选取2个元素是不可能的），所以整个表达式的结果为0。
因此， C ( 10 , 5 ) m o d 3 = 0 C(10, 5) \mod 3 = 0 C(10,5)mod3=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 a a， b b b,求 C a b C_{a}^{b} Cab的值。
注意结果可能很大，需要使用高精度计算。
输入格式
共一行，包含 a a a, b b b的值。
输出格式
共一行，输出 C a b C_{a}^{b} Cab的值
数据范围
1 ≤ b ≤ a ≤ 5000 1 \le b \le a \le 5000 1≤b≤a≤5000，
思路
这一题的区别是没有了模数，这意味着我们的答案会非常大，因此需要使用高精度来进行计算，我们直接通过计算式进行计算即可 C a b = a ! ( a − b ) ! b ! C_{a}^{b} = \frac{a!}{(a - b)!b!} Cab=(a−b)!b!a!，但是这样我们就需要进行高精度乘法以及高精度除法，因此可以进行一些优化。
优化的思路是进行质因数分解，需要统计出答案数的所有质因数的指数，这里涉及到的一个知识点就是一个数阶乘a!的某个质因子p的指数为 ⌊ a p ⌋ + ⌊ a p 2 ⌋ + ⌊ a p 3 ⌋ + ⋯ \lfloor \frac{a}{p} \rfloor + \lfloor \frac{a}{p^2} \rfloor + \lfloor \frac{a}{p^3} \rfloor + \cdots ⌊pa⌋+⌊p2a⌋+⌊p3a⌋+⋯。
这应该挺好理解的， a ! a! a!会涉及到所有的 1 ∼ a 1 \sim a 1∼a中的数，那么对于包含了 p p p作为质因子的这些数而言，我们首先统计最后能够贡献一次的数的个数，即 ⌊ a p ⌋ \lfloor \frac{a}{p} \rfloor ⌊pa⌋；但其中也会有能够贡献两次或更多次的数，也就是 p k ( k ≥ 2 ) p^k \; (k \ge 2) pk(k≥2)为其因子，对于这些数我们在统计一次 p 1 p^1 p1时已经统计了一次了，因此对于包含 p 2 p^2 p2为因子的数，只需再统计一次即可，即加上 ⌊ a p 2 ⌋ \lfloor \frac{a}{p^2} \rfloor ⌊p2a⌋，以此类推。
这里的代码会涉及到之前的筛质数以及高精度乘法，可以先去看一下。
代码
#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 n n个 0 0 0和 n n n个 1 1 1，它们将按照某种顺序排成长度为 2 n 2n 2n的序列，求他们能排列成的所有序列中，能构满足任意前缀序列中 0 0 0的个数都不少于 1 1 1的个数的序列有多少个。
输出的答案对 1 0 9 + 7 10^9+7 109+7取模
输入格式
共一行，包含整数 n n n。
输出格式
共一行，包含一个整数，表示答案。
数据范围
1 ≤ n ≤ 1 0 5 1 \le n \le 10^5 1≤n≤105，
思路
这一题其实是卡特兰数的一种数学模型，这里就先直接给出他的计算公式了（百度学了一下，应该不止这一种公式，这里就写y总讲的了）
h
(
n
)
=
C
2
n
n
−
C
2
n
n
−
1
=
C
2
n
n
n
+
1
h(n) = C_{2n}^{n} - C_{2n}^{n - 1} = \frac{C_{2n}^{n}}{n + 1}
h(n)=C2nn−C2nn−1=n+1C2nn
我们这里也用y总讲的模型进行讲解——可以将其看为一个路径计数的问题
我们将选
0
0
0看做向左走一步，选
1
1
1看做向上走一步，假设我们一共需要走12步，那么我们一半一半选的终点就是
(
6
,
6
)
(6,6)
(6,6)，并且所有的路径方案数为
C
2
n
n
=
C
12
6
C_{2n}^{n} = C_{12}^{6}
C2nn=C126下图中给出了一条符合题意的路径，大家可以看一下，到达终点前的所有前缀路径中向左走的步数都不小于向上走的步数，图中红线为对称线
要找出的是所有符合这种条件的路径，因此可以反向思考，也就是找到所有不符合条件的路径，在下图中给出了一条不符合题意的橙色路径，图中紫色斜线为另一条关键线，可以看到，当橙色路径第一次碰到紫色斜线时，表明其破坏了题目所要求的条件，我们也由此可以看出若一条路径会碰到或者穿过紫色斜线，其为不符合题意的路径，我们要做的就是找出所有这样的路径。
这里的做法比较巧妙，我们将橙色路径第一次碰到紫色斜线的部分后面进行轴对称，如下图绿色路径所示，所有经过紫色斜线的路径经过这样操作之后，其终点一定为点 ( 5 , 7 ) (5,7) (5,7)，因为在轴对称之前的终点固定是 ( 6 , 6 ) (6,6) (6,6)。
那么就可以得出这样的结论：所有从 ( 0 , 0 ) (0,0) (0,0)出发到达点 ( 5 , 7 ) (5,7) (5,7)的路径都可以由不满足题意的到达 ( 6 , 6 ) (6,6) (6,6)的路径经过上述变换得到，其数量为 C 2 n n − 1 = C 12 5 C_{2n}^{n - 1} = C_{12}^{5} C2nn−1=C125。
因此最后要求的路径数目就是
C
12
6
−
C
12
5
C_{12}^{6} - C_{12}^{5}
C126−C125。
最后实现到代码上我们肯定要进一步化简，也就是 h ( n ) = C 2 n n − C 2 n n − 1 = C 2 n n n + 1 h(n) = C_{2n}^{n} - C_{2n}^{n - 1} = \frac{C_{2n}^{n}}{n + 1} h(n)=C2nn−C2nn−1=n+1C2nn
看题目给出的数据范围是 1 0 5 10^5 105，因此需要使用的是上面组合数第二题的计算方法。
代码
#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;
}
![LosslessScaling-学习版[steam价值30元的游戏无损放大/补帧工具]](https://i-blog.csdnimg.cn/direct/a394c3998dca4e93b086110b9303607b.png)
