组合数学

1 计数原理和方法

1.1 加法原理

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

1.2 乘法原理

完成一件事情有 \(n\) 个步骤，第一个步骤有 \(n_1\) 个方法，第二个步骤有 \(n_2\) 个方法，\(\cdots\)，那么完成这件事共有 \(\prod\limits_{i=1}^nn_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)=\dfrac{n!}{(n-m)!}\)

1.3.2 组合

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

的组合，用 \(C_n^m\) 表示。

\(C_n^m=\dfrac{A_n^m}{A_m^m}=\dfrac{n!}{m!(n-m)!}\)

1.3.3 常用策略与方法

1.3.3.1 特殊优先

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

1.3.3.2 捆绑法

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

1.3.3.3 插空法

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

插空法常用结论：

将一个长度为 \(n\) 的序列划分为 \(m\) 非空序列的方案数为 \(C_{n-1}^{m-1}\)。
将一个长度为 \(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\) 个元素，有两种情况：

将它放在 \(n\) 位置，此时剩下的 \(n-1\) 个元素有 \(D(n-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=\dfrac{n!}{m!(n-m)!}\) 直接计算。

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

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

1.3.4.3 Lucas 定理

公式：

\[C_n^m\equiv C_{n/p}^{m/p}\times C_{n\bmod{p}}^{m\bmod{p}}\pmod{p} \]

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

时间复杂度为 \(O(p\log p+\log_pn)\)。

完整代码：

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_pn)\)	当 \(n\) 较大， \(p\) 为素数且 \(n>p\) 时使用

2 中国剩余定理（CRT）

2.1 简介

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

\[\begin{cases}x\equiv r_1\pmod{m_1}\\x\equiv r_2\pmod{m_2}\\\cdots\\x\equiv r_n\pmod{m_n}\end{cases} \]

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

2.2 算法流程

计算所有模数的乘积 \(M=\prod_{i=1}^nm_i\)。
对于第 \(i\) 个方程，计算出 \(c_i=\dfrac M{m_i}\)，并求出 \(c_i\) 在模 \(m_i\) 意义下的逆元。
方程的解为 \(x=\sum\limits_{i=1}^nr_ic_ic_i^{-1}\bmod{M}\)。

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

2.2.1 算法证明

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

先证明 \(\sum\limits_{i=1}^nr_ic_ic_i^{-1}\) 满足 \(x\equiv r_i\pmod{m_i}\)。

分情况讨论：

当 \(i\ne j\) 时，\(c_j\equiv0\pmod{m_i}\)，则 \(c_jc_j^{-1}\equiv0\pmod{m_i}\)（因为 \(c_j\) 中去掉了 \(m_j\) 但还有 \(m_i\)）

当 \(i=j\) 时，\(c_i\ne0\pmod{m_i}\)，则 \(c_ic_i^{-1}\equiv1\pmod{m_i}\)。

综上，对于任何一个 \(i\)，都有：

\(x\equiv\sum\limits_{j=1}^nr_jc_jc_j^{-1}\equiv r_ic_ic_i^{-1}\equiv r_i\pmod{m_i}\)

而 \(\bmod{M}\) 对于任何 \(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{cases}x\equiv r_1\pmod{m_1}\\x\equiv r_2\pmod{m_2}\\\cdots\\x\equiv r_n\pmod{m_n}\end{cases} \]

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

2.3.1 算法过程

考虑前两个方程 \(x\equiv r_1\pmod{m_1},x\equiv r_2\pmod{m_2}\)。

转化为不定方程即 \(x=m_1p+r_1=m_2q+r_2\)。

所以 \(m_1p-m_2q=r_2-r_1\)

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

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

于是通解就是 \(P = p+\dfrac {m_2}d*k,Q=q-\dfrac{m_1}d*k\)

所以 \(x=m_1P+r_1=\dfrac{m_1m_2}d*k+m_1p+r_1=\operatorname{lcm}(m_1m_2)*k+m_1p+r_1\)

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

其中 \(r=m_1p+r_1,m=\operatorname{lcm}(m_1m_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\pmod{p}\)，其中 \(p\) 不一定为素数。

3.2 求解

我们逐步分解问题求解。

3.2.1 CRT

将 \(p\) 进行质因数分解得 \(p=\prod\limits_{i}p_i^{k_i}\)。

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

3.2.2 组合数模素数幂

我们继续看上面的式子：\(C_n^m\bmod{p^k}\)。

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

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

\[\dfrac{\dfrac{n!}{p^{k_1}}}{\dfrac{m!}{p^{k_1}}\times\dfrac{(n-m)!}{p^{k_3}}}\times p^{k_1-k_2-k_3} \]

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

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

3.2.3 阶乘模素数幂

我们继续看下面的式子：\(\dfrac{n!}{p^a}\bmod{p^k}\)。

我们先算 \(n!\bmod{p^k}\)。

举个例子：\(n=22,p=3,k=2\)。

写出来：

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

把当中所有含 \(p\) 的数提出来，得：

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

可以看出上式分为三个部分：第一个部分是 \(3\) 的幂，次数是不大于 \(22\) 的 \(3\) 的倍数的个数，也就是 \(\lfloor\dfrac np\rfloor\)。

第二个部分是 \(7!\)，也就是 \(\lfloor\dfrac np\rfloor!\)，递归求解。

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

\(1*2*4*5*7*8\equiv10*11*13*14*16*17\pmod{p^k}\)。

写成如下形式更加显然。

\[\prod_{i,\gcd(i,p)=1}^{p^k}i\equiv\prod_{i,\gcd(i,p)=1}^{p^k}(i+tp^k)\pmod{p^k}(t\in\mathbb{N^*}) \]

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

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

然后三部分的乘积就是 \(n!\)，我们继续看 \(n\bmod{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 的式子