组合数学

合集 - 学习笔记(52)

1.随机化算法2024-02-27

2.组合数学2024-02-27

3.概率和期望2024-02-274.树链剖分2024-02-275.矩阵乘法与矩阵快速幂2024-02-276.最短路2024-02-277.卡特兰数、Prüfer 序列、BSGS2024-02-278.分块和莫队2024-02-279.AC 自动机2024-02-2710.平衡树2024-02-2711.基础字符串算法2024-02-2712.基础数论2024-02-2713.Miller-Rabin 与 Pollard-Rho2024-08-2314.广义后缀自动机2024-06-2515.后缀自动机 SAM2024-06-2216.回文自动机 PAM2024-06-1817.Manacher2024-06-1618.后缀数组2024-06-0819.01 分数规划2024-06-0220.网络流2024-05-0521.整体二分2024-05-0522.cdq 分治2024-05-0423.点分治2024-04-2124.虚树2024-04-1425.扫描线2024-04-1126.笛卡尔树2024-04-0627.基环树2024-04-0228.树哈希2024-03-0929.莫比乌斯反演2024-11-2830.二分图2024-11-0731.朱刘算法2024-10-2932.最大团2024-10-2933.杂项2024-08-2434.杜教筛2024-12-0135.拉格朗日插值2024-12-1036.线段树综合2024-12-1737.可持久化数据结构2024-12-2438.K-D Tree2024-12-2539.Burnside 引理与 Polya 定理2024-12-2740.线性基2024-12-3041.替罪羊树01-0942.LCT01-1843.插头 dp01-2044.原根01-2045.多项式乘法01-2246.斯特林数01-2347.二项式反演与斯特林反演01-2348.Min-Max 容斥01-2449.辛普森积分法01-2550.Min_25 筛02-0351.凸包02-1652.2-SAT03-03

收起

1 计数原理和方法

1.1 加法原理

完成一件事情有 $n$ 个办法，第一类方法有 $n_1$ 个方案，第二类方法有 $n_2$ 个方案，$\cdots$，那么完成这件事共有 $\sum _{i=1}^n{n}_i$ 种方法。

1.2 乘法原理

完成一件事情有 $n$ 个步骤，第一个步骤有 $n_1$ 个方法，第二个步骤有 $n_2$ 个方法，$\cdots$，那么完成这件事共有 $\prod _{i=1}^n{n}_i$ 种方法。

1.3 排列与组合

1.3.1 排列

从 $n$ 个不同的元素中选出 $m$ 个，按照一定顺序排成一列，叫做从 $n$ 个不同的元素中选出 $m$ 个

的排列，用 $A_n^m$ 表示。

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

1.3.2 组合

从 $n$ 个不同的元素中选出 $m$ 个，不考虑顺序，叫做从 $n$ 个不同的元素中选出 $m$ 个

的组合，用 $C_n^m$ 表示。

$C_n^m={\displaystyle \frac{A_n^m}{A_m^m}}={\displaystyle \frac{n!}{m!(n-m)!}}$

1.3.3 常用策略与方法

1.3.3.1 特殊优先

将特殊要求的优先考虑，普通情况靠后考虑。

1.3.3.2 捆绑法

将要求相邻的元素捆绑为一个整体，与其他元素算出排列，再在整体内部进行排列。

1.3.3.3 插空法

插空法用于处理不相邻问题。把要求不相邻的元素放一边，算出其他元素的排列，再把不相邻的元素插入已经排好的元素之间的空位。

插空法常用结论：

1. 将一个长度为 $n$ 的序列划分为 $m$ 非空序列的方案数为 $C_{n-1}^{m-1}$。
2. 将一个长度为 $n$ 的序列划分为 $m$ 可空序列的方案数为 $C_{m+n-1}^{m-1}$。

1.3.3.4 倍缩法

对于某几个元素顺序一定的排列，先把这几个元素和其他的一起排列，再除以这几个元素间的全排列数。

1.3.3.5 环问题策略

把环破为链，而由于圆没有首尾之分。所以例如 $123456$ 的排列与 $234561$ 和 $345612$ 本质上是相同的，所以还要除以 $n$。

一般的，$n$ 个元素做环排列，有 $(n-1)!$ 种方法。

1.3.3.6 错排问题

错排问题指的是有 $n$ 个元素，若一个序列的所有数都不在他原来的位置上，就说这是原序列的一个错排。

这是一个递推问题。将 $n$ 个数的错排个数记作 $D(n)$，然后分步分类讨论。

第一步，考虑放下第 $n$ 个元素，则有 $n-1$ 种方法。设当前放下的位置是 $k$。

第二步，考虑放下第 $k$ 个元素，有两种情况：

1. 将它放在 $n$ 位置，此时剩下的 $n-1$ 个元素有 $D(n-2)$ 种方法。
2. 将它不放在 $n$ 位置，这时对于除 $n$ 外的 $n-1$ 个元素有 $D(n-1)$ 种方法。

综上，$D(n)=(n-1)\times (D(n-2)+D(n-1))$。

边界为 $D(1)=0,D(2)=1$。

1.3.4 组合数取模

1.3.4.1 递推求解

递推式明显可得：$C_n^m=C_{n-1}^m+C_{n-1}^{m-1}$。加上取模即可。

时间复杂度 $O(n^2)$，太劣。

1.3.4.2 公式法求解

利用 $C_n^m={\displaystyle \frac{n!}{m!(n-m)!}}$ 直接计算。

用两个数组存储模 $p$ 意义下阶乘和阶乘的逆元即可计算。

补：阶乘的逆元可以递推计算 $(n!{)}^{-1}\equiv i^{p-2}\times (n-1){!}^{-1}(\mathrm{mod}p)$。

1.3.4.3 Lucas 定理

公式：

$$C_n^m\equiv C_{n/p}^{m/p}\times C_{nmodp}^{mmodp}(\mathrm{mod}p)$$

其中第一项继续用 Lucas 定理递归求解，第二项直接用公式法求解即可。

时间复杂度为 $O(p\log p+{\log }_{p}n)$。

完整代码：

int g[Maxn], f[Maxn], p;

int qpow(int a, int b) {
	int res = 1;
	while(b) {
		if(b & 1) {
			res = (res * a) % p;
		}
		a = (a * a) % p;
		b >>= 1;
	}
	return res;
}

void init() {
	f[0] = g[0] = 1;
	for(int i = 1; i < p; i++) {
		f[i] = (f[i - 1] * i) % p;
		g[i] = (g[i - 1] * qpow(i, p - 2)) % p;
	}
}

int getc(int n, int m) {
	if(n < m) return 0;
	return f[n] * g[n - m] % p * g[m] % p;
}

int lucas(int n, int m) {
	if(m == 0) return 1;
	if(n < m) return 0;
	return lucas(n / p, m / p) * getc(n % p, m % p) % p;
}

1.3.4.4 小结

比较三种方法：

方法	复杂度	用途
递推法	$O(n^2)$	当 $n$ 较小时可以使用
公式法	预处理 $O(p)$，单次询问 $O(1)$	当 $n$ 较大，$p$ 为素数且 $n<p$ 时使用
Lucas 定理	$O(p\log p+{\log }_{p}n)$	当 $n$ 较大， $p$ 为素数且 $n>p$ 时使用

2 中国剩余定理（CRT）

2.1 简介

中国剩余定理，简称 CRT，用来解如下的一元同余方程组的最小非负数解 $x$：

$$\{\begin{array}{l}x\equiv r_1(\mathrm{mod}{m}_1)\\ x\equiv r_2(\mathrm{mod}{m}_2)\\ \cdots \\ x\equiv r_n(\mathrm{mod}{m}_n)\end{array}$$

（其中 $m_1,m_2,\cdots ,m_n$ 两两互质）

2.2 算法流程

1. 计算所有模数的乘积 $M=\prod _{i=1}^n{m}_i$。
2. 对于第 $i$ 个方程，计算出 $c_i={\displaystyle \frac{M}{m_i}}$，并求出 $c_i$ 在模 $m_i$ 意义下的逆元。
3. 方程的解为 $x=\sum _{i=1}^n{r}_i{c}_i{c}_i^{-1}modM$。

复杂度 $O(n\log c)$。

2.2.1 算法证明

虽然没用，但是理解了更好记忆吧。

先证明 $\sum _{i=1}^n{r}_i{c}_i{c}_i^{-1}$ 满足 $x\equiv r_i(\mathrm{mod}{m}_i)$。

分情况讨论：

当 $i\ne j$ 时，$c_j\equiv 0(\mathrm{mod}{m}_i)$，则 $c_j{c}_j^{-1}\equiv 0(\mathrm{mod}{m}_i)$（因为 $c_j$ 中去掉了 $m_j$ 但还有 $m_i$）

当 $i=j$ 时，$c_i\ne 0(\mathrm{mod}{m}_i)$， 则 $c_i{c}_i^{-1}\equiv 1(\mathrm{mod}{m}_i)$。

综上，对于任何一个 $i$，都有：

$x\equiv \sum _{j=1}^n{r}_j{c}_j{c}_j^{-1}\equiv r_i{c}_i{c}_i^{-1}\equiv r_i(\mathrm{mod}{m}_i)$

而 $modM$ 对于任何 $m_i$ 来说只是减去了 $m_i$ 的若干倍，对余数 $r_i$ 没有影响。

算法证毕。

2.2.2 代码

直接模拟即可。

int exgcd(int a, int b, int &x, int &y) {//扩欧求逆元
	if(b == 0) {
		x = 1, y = 0;
		return a;
	}
	int ret = exgcd(b, a % b, x, y);
	int t = x;
	x = y;
	y = t - a / b * y;
	return ret;
}

int CRT(int m[], int r[]) {
	int M = 1, ans = 0;
	for(int i = 1; i <= n; i++) {
		M *= m[i];
	}
	for(int i = 1; i <= n; i++) {
		int c = M / m[i], x, y;
		int d = exgcd(c, m[i], x, y);
		ans = (ans + r[i] * c * x % M) % M;
	}
	ans = (ans % M + M) % M;
	return ans;
}

2.3 扩展中国剩余定理（exCRT）

扩展中国剩余定理，简称 exCRT，用来解如下的一元同余方程组的最小非负数解 $x$：

$$\{\begin{array}{l}x\equiv r_1(\mathrm{mod}{m}_1)\\ x\equiv r_2(\mathrm{mod}{m}_2)\\ \cdots \\ x\equiv r_n(\mathrm{mod}{m}_n)\end{array}$$

其中 $m_1,m_2,\cdots ,m_n$ 不一定两两互质。

2.3.1 算法过程

考虑前两个方程 $x\equiv r_1(\mathrm{mod}{m}_1),x\equiv r_2(\mathrm{mod}{m}_2)$。

转化为不定方程即 $x=m_{1}p+r_1=m_{2}q+r_2$。

所以 $m_{1}p-m_{2}q=r_2-r_1$

由裴蜀定理可知当 $gcd(m_1,m_2)\mid (r_2-r_1)$ 时有解。

由扩欧，得到方程特解为 $p=p\ast {\displaystyle \frac{r_2-r_1}{d}},q=q\ast {\displaystyle \frac{r_2-r_1}{d}}$（$d$ 即为 $gcd(m_1,m_2)$）

于是通解就是 $P=p+{\displaystyle \frac{m_2}{d}}\ast k,Q=q-{\displaystyle \frac{m_1}{d}}\ast k$

所以 $x=m_{1}P+r_1={\displaystyle \frac{m_1{m}_2}{d}}\ast k+m_{1}p+r_1=\mathrm{lcm}(m_1{m}_2)\ast k+m_{1}p+r_1$

根据上式，前两个方程就等价于合并为一个方程：$x\equiv r(\mathrm{mod}m)$

其中 $r=m_{1}p+r_1,m=\mathrm{lcm}(m_1{m}_2)$。

于是 $n$ 个方程合并 $n-1$ 次即可求解。

2.3.2 代码

int exgcd(int a, int b, int &x, int &y) {
	if(b == 0) {
		x = 1, y = 0;
		return a;
	}
	int ret = exgcd(b, a % b, x, y);
	int t = x;
	x = y;
	y = t - a / b * y;
	return ret;
}

int exCRT(int m[], int r[]) {
	int m1, m2, r1, r2, p, q;
	m1 = m[1], r1 = r[1];
	for(int i = 2; i <= n; i++) {
		m2 = m[i], r2 = r[i];
		int d = exgcd(m1, m2, p, q);
		if((r2 - r1) % d != 0) return -1;
		p = p * (r2 - r1) / d;
		p = (p % (m2 / d) + (m2 / d)) % (m2 / d);
		r1 = m1 * p + r1;
		m1 = m1 * m2 / d;
	}
	return (r1 % m1 + m1) % m1;
}

3 扩展卢卡斯定理（exLucas）

虽然这应该是放在 1.3.4 中的，但是他值得单独拿出来。

我愿称之为数论与组合的集大成之算法，他使用到的包括 CRT、exgcd、快速幂、分解质因数。

3.1 问题

求 $C_n^m(\mathrm{mod}p)$，其中 $p$ 不一定为素数。

3.2 求解

我们逐步分解问题求解。

3.2.1 CRT

将 $p$ 进行质因数分解得 $p=\prod _i{p}_i^{k_i}$。

此时显然 $p_i$ 两两互质，所以如果求出所有 $C_n^{m}mod{p}_i^{k_i}=a_i$，那么就可以构造出若干个类似 $C_n^m=a_i(\mathrm{mod}{p}_i^{k_i})$，此时利用中国剩余定理即可求解出结果。

3.2.2 组合数模素数幂

我们继续看上面的式子：$C_n^{m}mod{p}^k$。

我们回顾组合数的公式 $C_n^m={\displaystyle \frac{n!}{m!(n-m)!}}$，发现由于 $m!$ 和 $(n-m)!$ 中可能有质因子 $p$，所以无法直接求逆元。

我们可以将 $n!,m!,(n-m)!$ 中的所有质因子 $p$ 都提出来，最后乘回去即可。也就可以变为下面式子：

$${\displaystyle \frac{{\displaystyle \frac{n!}{p^{k_1}}}}{{\displaystyle \frac{m!}{p^{k_1}}}\times {\displaystyle \frac{(n-m)!}{p^{k_3}}}}}\times p^{k_1-k_2-k_3}$$

此时分母就互质了，直接求逆元即可。

但是还有一个问题：如何求阶乘模数。

3.2.3 阶乘模素数幂

我们继续看下面的式子：${\displaystyle \frac{n!}{p^a}}mod{p}^k$。

我们先算 $n!mod{p}^k$。

举个例子：$n=22,p=3,k=2$。

写出来：

$22!=1\ast 2\ast 3\ast 4\ast 5\ast 6\ast 7\ast 8\ast 9\ast 10\ast 11\ast 12\ast 13\ast 14\ast 15\ast 16\ast 17\ast 18\ast 19\ast 20\ast 21\ast 22$

把当中所有含 $p$ 的数提出来，得：

$22!=3^7\ast (1\ast 2\ast 3\ast 4\ast 5\ast 6\ast 7)\ast (1\ast 2\ast 4\ast 5\ast 7\ast 8\ast 10\ast 11\ast 13\ast 14\ast 16\ast 17\ast 19\ast 20\ast 22)$

可以看出上式分为三个部分：第一个部分是 $3$ 的幂，次数是不大于 $22$ 的 $3$ 的倍数的个数，也就是 $\lfloor {\displaystyle \frac{n}{p}}\rfloor$。

第二个部分是 $7!$，也就是 $\lfloor {\displaystyle \frac{n}{p}}\rfloor !$，递归求解。

第三个部分是 $n!$ 中与 $p$ 互质的部分之积，这一部分具有如下性质：

$1\ast 2\ast 4\ast 5\ast 7\ast 8\equiv 10\ast 11\ast 13\ast 14\ast 16\ast 17(\mathrm{mod}{p}^k)$。

写成如下形式更加显然。

$$\prod _{i,gcd(i,p)=1}^{p^k}i\equiv \prod _{i,gcd(i,p)=1}^{p^k}(i+t{p}^k)(\mathrm{mod}{p}^k)(t\in {\mathbb{N}}^{\mathbb{\ast }})$$

暴力求出 $\prod _{i,gcd(i,p)=1}^{p^k}i$ 然后用快速幂求它的 $\lfloor {\displaystyle \frac{n}{p^k}}\rfloor$ 次幂。

最后再乘上 $19\ast 20\ast 22$，一般形式为 $\prod _{i,gcd(i,p)=1}^{nmod{p}^k}i$，暴力乘上去即可。

然后三部分的乘积就是 $n!$，我们继续看 $nmod{p}^a$，此时第一部分全部消掉了，而第二部分递归计算时已经除掉了质因子 $p$。因此，最后结果为第二部分的递归结果与第三部分的乘积。

直接看代码：

int fac(int n, int p, int mod) {//阶乘模素数幂
	if(!n) return 1;
	int ans = 1;
	for(int i = 1; i < mod; i++) { 
		if(i % p) {
			ans = ans * i % mod;
		}
	}
	ans = qpow(ans, n / mod, mod);
	for(int i = 1; i <= n % mod; i++) {
		if(i % p) {
			ans = ans * i % mod;
		}
	}//全部是第三部分
	return ans * fac(n / p, p, mod)/*递归求解第二部分*/ % mod; 
}

3.2.4 回代求解——组合数模素数幂

回到 3.2.2 的式子

$${\displaystyle \frac{{\displaystyle \frac{n!}{p^{k_1}}}}{{\displaystyle \frac{m!}{p^{k_1}}}\times {\displaystyle \frac{(n-m)!}{p^{k_3}}}}}\times p^{k_1-k_2-k_3}$$

现在就可以直接转化为代码了：

int C(int n, int m, int p, int mod) {//组合数模素数幂
	if(n < m) {
		return 0;
	}
	int f1 = fac(n, p, mod), f2 = fac(m, p, mod), f3 = fac(n - m, p, mod);
	//求出除掉 p 后的阶乘值 
	int cnt = 0;//k1-k2-k3
	for(int i = n; i; i /= p) {
		cnt += i / p;
	}
	for(int i = m; i;i i /= p) {
		cnt -= i / p;
	}
	for(int i = n - m; i; i /= p) {
		cnt -= i / p;
	}
	return f1 * inv(f2, mod) % mod * inv(f3, mod) % mod * qpow(p, cnt, mod) % mod;
	//利用逆元计算 
} 

3.2.5 回代求解——CRT

直接带入计算即可。

代码直接看模板吧。

3.3 代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

typedef long long LL;
const int Maxn = 2e6 + 5;

int qpow(int a, int b, int mod) {//快速幂 
	int res = 1;
	while(b) {
		if(b & 1) {
			res = (res * a) % mod;
		}
		a = (a * a) % mod;
		b >>= 1;
	}
	return res;
}

int exgcd(int a, int b, int &x, int &y) {//扩欧 
	if(!b) {
		x = 1, y = 0;
		return a;
	}
	int ret = exgcd(b, a % b, x, y);
	int t = x;
	x = y;
	y = t - a / b * y;
	return ret;
}

int inv(int a, int p) {//利用扩欧求逆元 
	int x, y;
	int ret = exgcd(a, p, x, y);
	return (x % p + p) % p;
}

int fac(int n, int p, int mod) {//阶乘模素数幂
	if(!n) return 1;
	int ans = 1;
	for(int i = 1; i < mod; i++) { 
		if(i % p) {
			ans = ans * i % mod;
		}
	}
	ans = qpow(ans, n / mod, mod);
	for(int i = 1; i <= n % mod; i++) {
		if(i % p) {
			ans = ans * i % mod;
		}
	}//全部是第三部分
	return ans * fac(n / p, p, mod)/*递归求解第二部分*/ % mod; 
}

int C(int n, int m, int p, int mod) {//组合数模素数幂
	if(n < m) {
		return 0;
	}
	int f1 = fac(n, p, mod), f2 = fac(m, p, mod), f3 = fac(n - m, p, mod);
	//求出除掉 p 后的阶乘值 
	int cnt = 0;//k1-k2-k3
	for(int i = n; i; i /= p) {
		cnt += i / p;
	}
	for(int i = m; i; i /= p) {
		cnt -= i / p;
	}
	for(int i = n - m; i; i /= p) {
		cnt -= i / p;
	}
	return f1 * inv(f2, mod) % mod * inv(f3, mod) % mod * qpow(p, cnt, mod) % mod;
	//利用逆元计算 
} 

int CRT(int m[], int n[], int l) {//普通 CRT
	int M = 1, ans = 0;
	for(int i = 1; i <= l; i++) {
		M *= m[i];
	}
	for(int i = 1; i <= l; i++) {
		int t = M / m[i];
		int p = inv(t, m[i]);
		ans = (ans + n[i] * t % M * p % M) % M;
	}
	ans = (ans % M + M) % M;
	return ans;
}

int a[Maxn], b[Maxn], cnt;

int exLucas(int n, int m, int p) {
	int w = sqrt(p);
	for(int i = 2; i <= w; i++) {
		int tmp = 1;
		while(p % i == 0) {//直接求 p^k 
			p /= i;
			tmp *= i;
		}
		if(tmp > 1) {
			a[++cnt] = C(n, m, i, tmp);//求 Cnm mod p^k 
			b[cnt] = tmp;
		}
	}
	if(p > 1) {
		a[++cnt] = C(n, m, p, p);
		b[cnt] = p;
	}
	return CRT(b, a, cnt);//CRT 求解 
}

int n, m, p;

signed main() {
	ios::sync_with_stdio(0);
	cin >> n >> m >> p;
	cout << exLucas(n, m, p);
	return 0;
}

长达 119 行，甚至超过了普通线段树。