组合数学复习笔记

· · 算法·理论

加法 & 乘法原理

这两者的区别是:加法分类,乘法分步。

排列组合

排列

n不同元素中取 m 个元素排成有序的一列,方案数用 A_{n}^{m} 表示。

我们这样计算:第一个数有 n 种选法,第二个数有 n-1 种选法……第 m 个数有 n-m+1 种选法,由乘法原理得

A_{n}^{m}=n(n-1)(n-2) \cdots (n-m+1)=\prod_{i=n-m+1}^{n} i=\dfrac{n!}{(n-m)!}

规定 m>nA_{n}^{m}=0

而如果 m=n,也就是所有元素都参与排列,此时称为全排列。全排列的计算

A_{n}^{n}=n!

全排列的枚举

在 C++ 标准库中提供了 next_permutation 函数,可以这样使用:

// a为要枚举排列的数组,n为长度
sort(a + 1, a + n + 1);
do {
  ...
} while (next_permutation(a + 1, a + n + 1));

时间复杂度为 O(n!)

有重复元素的排列问题

顾名思义就是元素有重,此时需要除去重复的排列。

考虑一组有 m 个的重复元素,其造成重复的个数就是它在排列中的次序交换,也就是 m! 种情况,所以总排列数

\dfrac{n!}{m_1!m_2!m_3! \cdots}

其中 n 为总元素个数,m 为各元素出现次数。

圆排列

n不同元素中取 m 个元素排成有序的一圈,方案数用 Q_{n}^{m} 表示。考虑断环为链,共有 m 个断点,故

Q_{n}^{m}=\dfrac{A_{n}^{m}}{m}=\dfrac{n!}{m(n-m)!}

组合

n不同元素中取 m 个元素组成无序的一个集合,方案数用 C_{n}^{m} 表示。

考虑先作排列,由于集合无序,需要除去重复的组合,也就是 m 的全排列,逆用乘法原理:

C_{n}^{m}=\dfrac{A_{n}^{m}}{m!}=\dfrac{n!}{m!(n-m)!}

规定 m>nC_{n}^{m}=0

组合的枚举

这是 pj 组题目,代码如下:

int n, m, cur[N];
void dfs(int x) {
    if (x > m) return ..., void();  // 这里cur就是一个组合,进行处理
    for (int i = cur[x - 1] + 1; i <= n; i++) cur[x] = i, dfs(x + 1);
}

dfs(1);

它用于枚举 1n 所有自然数选 m 个的组合。如果套上 next_permutation,还能实现排列的枚举。

二项式定理

(a+b)^n=\sum_{i=0}^{n} C_{n}^{i} a^{n-i} b^{i}

这个定理表明,二项式 (a+b)^n 展开项的系数与组合数有直接关系。

我们在初中阶段学习过杨辉三角:

其每一个位置的数字通过左上 + 右上确定,每一行所对应的就是二项式展开的系数。

组合数的性质

  1. 对称性:
C_{n}^{m} = C_{n}^{n-m}

代数推导易证,组合意义就是把选出的集合取补集。

  1. 递推式:
C_{n}^{m}=C_{n-1}^{m}+C_{n-1}^{m-1}

代数推导略去,组合意义类似 dp 的思想可以证明。这个式子实际上就是杨辉三角的递推式。

  1. 二项式定理的特殊情况:
2^n=\sum_{i=0}^{n} C_n^i

也就是杨辉三角每一行的和。

  1. 斐波那契数列:
fib_{n+1}=\sum_{i=0}^{n} C_{n-i}^i

把杨辉三角每条斜 30 \degree 角的线取出来相加可以发现。

  1. 范德蒙恒等式:
C_{m+n}^{k}=\sum_{i=0}^{k} C_m^i C_n^{k-i}

假设有两堆物品,每堆分别有 m,n 个物品,总共取 k 个,则方案数可分解为:从第一堆取 i 个物品,第二堆取 k-i 个物品,且两种选择独立,故由乘法得到贡献。最后方案即为求和。

计算组合数

定义法

根据定义直接计算,注意先乘后除。

inline long long C(int n, int m) {
    if (m > n) return 0;
    long long res = 1;
    for (int i = 1; i <= m; i++) {
        res = res * (n - i + 1) % mod * power(i, mod - 2) % mod;
    }
    return res;
}

递推法

即利用组合数性质 2 预处理杨辉三角。

for (int i = 0; i < N; i++) c[i][0] = 1, c[i][i] = 1;
for (int i = 0; i < N; i++) {
  for (int j = 1; j < i; j++) c[i][j] = c[i - 1][j] + c[i - 1][j - 1];
}

逆元法

通过预处理阶乘和阶乘的逆元,直接用定义式计算组合数。

inline long long power(long long a, long long b) {
    long long res = 1;
    for (a %= mod; b; b >>= 1) {
        if (b & 1) res = res * a % mod;
        a = a * a % mod;
    }
    return res;
}

long long fac[N], ifac[N];
inline void pre() {
    fac[0] = fac[1] = ifac[0] = ifac[1] = 1;
    for (int i = 2; i < N; i++) fac[i] = fac[i - 1] * i % mod, ifac[i] = power(fac[i], mod - 2);
}

inline long long C(int n, int m) {
    if (n < m) return 0;
    return fac[n] * ifac[m] % mod * ifac[n - m] % mod; 
}

Lucas 定理法

Lucas 定理如下:

C_n^m \equiv C_{n \bmod p}^{m \bmod p} C_{\lfloor \frac{n}{p}\rfloor}^{\lfloor \frac{m}{p}\rfloor} \pmod p

其中 p 是质数。利用这个式子,先用逆元法预处理所有 p 以内的组合数,查询时递归计算即可。

inline long long power(long long a, long long b) {
    long long res = 1;
    for (a %= mod; b; b >>= 1) {
        if (b & 1) res = res * a % mod;
        a = a * a % mod;
    }
    return res;
}

long long fac[N], ifac[N];
inline void pre() {
    fac[0] = fac[1] = ifac[0] = ifac[1] = 1;
    for (int i = 2; i < N; i++) fac[i] = fac[i - 1] * i % mod, ifac[i] = power(fac[i], mod - 2);
}

inline long long C(int n, int m) {
    if (n < m) return 0;
    return fac[n] * ifac[m] % mod * ifac[n - m] % mod; 
}
long long lucas(long long n, long long m) {
    if (m == 0) return 1;
    if (n < mod && m < mod) return C(n, m);
    return lucas(n / mod, m / mod) * C(n % mod, m % mod) % mod;
}

P3807 【模板】卢卡斯定理/Lucas 定理。

exLucas 法

exLucas 跟 Lucas 没有半毛钱关系。

将模数 p 质因数分解为 p=p_1^{a_1} p_2^{a_2} \cdots,然后计算 C_n^m \bmod p^a,最后再把答案用 exCRT 合并。

现在的重点是算 C_n^m \bmod p^a,由定义式得所求即为 \dfrac{n!}{m!(n-m)!} \bmod p^a,只要求出阶乘及逆元即可。先将阶乘中 p 的倍数算掉,即 \lfloor \dfrac{n}{p}\rfloor!,而剩余项存在循环节,可以一起计算,凑不进循环节的单独算。

For example:

\begin{aligned} 20!&=1\times 2\times 3\times 4\times 5\times 6\times 7\times 8\times 9\times 10\times 11\times 12\times 13\times 14\times 15\times 16\times 17\times 18\times 19\times 20\\ &=(1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19\times 20) \times 6! \times 3^6\\ &=(1\times 2 \times 4 \times 5 \times 7 \times 8)^2\times 19 \times 20 \times 6! \times 3^6 \end{aligned}

至于阶乘逆元,用 exGCD 求得即可。

inline long long power(long long a, long long b, long long p) {
    long long res = 1;
    for (a %= p; b; b >>= 1) {
        if (b & 1) res = res * a % p;
        a = a * a % p;
    }
    return res;
}

void exgcd(long long a, long long b, long long& x, long long& y) {
    if (b == 0) return x = 1, y = 0, void();
    exgcd(b, a % b, x, y);
    int t = x;
    x = y, y = t - a / b * y;
}
inline long long inverse(long long x, long long p) {
    long long a, b;
    exgcd(x, p, a, b);
    return (a % p + p) % p;
}

inline long long calc(long long n, long long x, long long p) {
    if (n == 0) return 1;
    long long s = 1;
    for (long long i = 1; i <= p; i++) {
        if (i % x) s = s * i % p;
    }
    s = power(s, n / p, p);
    for (long long i = n / p * p + 1; i <= n; i++) {
        if (i % x) s = i % p * s % p;
    }
    return s * calc(n / x, x, p) % p;
}

inline long long multilucas(long long m, long long n, long long x, long long p) {
    int cnt = 0;
    for (long long i = m; i; i /= x) cnt += i / x;
    for (long long i = n; i; i /= x) cnt -= i / x;
    for (long long i = m - n; i; i /= x) cnt -= i / x;
    return power(x, cnt, p) % p * calc(m, x, p) % p * inverse(calc(n, x, p), p) % p
           * inverse(calc(m - n, x, p), p) % p;
}

long long x[20];
inline int crt(int n, long long* a, long long* m) {
    long long mod = 1, res = 0;
    for (int i = 1; i <= n; i++) mod *= m[i];
    for (int i = 1; i <= n; i++) {
        x[i] = mod / m[i];
        long long x0 = -1, y0 = -1;
        exgcd(x[i], m[i], x0, y0);
        res += a[i] * x[i] * (x0 >= 0 ? x0 : x0 + m[i]);
    }
    return res % mod;
}

inline long long exlucas(long long m, long long n, long long p) {
    int len = 0;
    long long p0[20], a0[20];
    for (long long i = 2; i * i <= p; i++) {
        if (p % i) continue;
        p0[++len] = 1;
        while (p % i == 0) p0[len] *= i, p /= i;
        a0[len] = multilucas(m, n, i, p0[len]);
    }
    if (p > 1) p0[++len] = p, a0[len] = multilucas(m, n, p, p);
    return crt(len, a0, p0);
}

P4720 【模板】扩展卢卡斯定理/exLucas。

经典模型

捆绑法

Q: 有 n+m 个不同元素要进行排列,其中 m 个元素必须连续,求方案数。

A: 将 m 个元素捆绑在一起,内部排列方案为 m! 种,这 m 个元素的整体视为一个元素,则外部排列方案为 (n+1)! 种,由乘法原理得方案数为 (n+1)!m!

插空法

Q: 有 n+m 个不同元素要进行排列,其中 m 个元素必须两两不相邻,求方案数。

A: 先将 n 个元素全排列,再求 n+1 个空隙而插入 m 块板子的方案数为 A_{n+1}^{m},所以方案数为 n! A_{n+1}^{m}

插板法

Q1: 现有 n 个相同元素,分为 k 组,每组至少有一个元素,求方案数。
(可以抽象为:求方程 x_1+x_2+\cdots+x_k=n整数解数目)

A1: 不考虑分组,而是考虑间断点,将 k-1 块板子插入到 n-1 个空里,则答案为 C_{n-1}^{k-1}

Q2: 现有 n 个相同元素,分为 k 组,每组可以为空,求方案数。
(可以抽象为:求方程 x_1+x_2+\cdots+x_k=n非负整数解数目)

A2: 先借 k 个元素过来放到每组里,由 Q1 可得方案数为 C_{n+k-1}^{k-1}=C_{n+k-1}^{n},再把 k 个元素拿走,方案数不变。

Q3: 现有 n 个相同元素,分为 k 组,第 i 组的元素数目不小于 a_i,求方案数。
(可以抽象为:求方程 x_1+x_2+\cdots+x_k=n 的解数目,其中 x_i \ge a_i

A3: 先借 \sum a_i 个元素过来,令 x_i'=x_i-a_i,可知 \sum x_i'=n-\sum a_i,由 Q2 得答案为 C_{n-\sum a_i+k-1}^{n-\sum a_i}

错位排列

错排数 D_i 是指将 1n 的自然数排列重新排列为 P_i 后,满足 \forall i \in [1,n], P_i \ne i 的方案数。

例如,n=3 时,错位排列有 \{2,3,1\}\{3,1,2\}

计算错排数可以递推。考虑 D_n 时,先把 n 放在 P_n,然后有两种情况:

  1. 前面 n-1 个数已经是错位排列;
  2. 前面 n-1 个数有一个在原位上,其他错位。

对于情况 1,第 n 个数可以与任一数字交换,有 (n-1)D_{n-1} 种方案;

对于情况 2,第 n 个数只能与原位上的交换,有 (n-1)D_{n-2} 种方案;

因此,

D_n=(n-1)(D_{n-1}+D_{n - 2})

另外地,D_1=0,D_2=1

long long D[N];
inline void calc() {
    D[1] = 0, D[2] = 1;
    for (int i = 3; i < N; i++) D[i] = (i - 1) * (D[i - 1] + D[i - 2]);
}

板子题为 P1595 信封问题。

Catalan 数

Catalan 数的起源是凸多边形三角剖分问题,即对于一个凸 n 边形,有多少种方式可以用不相交的对角线将其划分为若干个三角形。这里将方案数记作 H_n。Catalan 数列为:

1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862\cdots

解决问题

  1. 括号序列计数

Q: 计算 n 对括号能形成多少种合法括号序列。

A: 方案数为 H_n。因为序列的第一个括号必须为 (,设之后有一个长度为 m 的子序列,然后为 ),再然后是另一个子序列,长度为 n-m-1,则方案为 H_n=\sum_{m=0}^{n-1} H_i H_{n-m-1},这就是 Catalan 的递推式。

  1. 不同构二叉树计数

Q: 计算 n 个节点的不同构二叉树数目。

A: 数目为 H_n。设根节点左子树大小为 m,则右子树大小为 n-m-1,左右子树又是一个子问题,所以方案为 H_n=\sum_{m=0}^{n-1} H_i H_{n-m-1},符合递推式。

  1. 三角剖分问题

Q: 求对于一个凸 n 边形,有多少种方式可以用不相交的对角线将其划分为若干个三角形。

A: 方案数为 H_n。固定多边形的一条边并选择一个与它不共线的初始顶点构成一个三角形,这个三角形将原多边形分为两个小多边形,大小分别为 mn-m-1,所以方案还是那个 H_n=\sum_{m=0}^{n-1} H_i H_{n-m-1}

  1. 不越过对角线的网格路径

Q: 在 n \times n 网格中,从 (0,0) 走到 (n,n) 且只能向右或向上移动,求不越过对角线的路径数。

A: 路径数为 H_n。证明可见 this。

  1. 出栈序列计数

Q: 入栈序列为一个 1n 的排列,求合法出栈序列数目。

A: 数目为 H_n。设第一个入栈元素在第 m 个出栈,则前面 m-1 个必须在之前出入栈,而之后又 n-m-1 种,递推式又双叒叕是 H_n=\sum_{m=0}^{n-1} H_i H_{n-m-1}

由此可见,能用 Catalan 数解决的问题都满足一个递归分解结构:一个大小为 n 的问题可以分解为两个独立子问题,规模分别为 mn-m-1,且两子问题分步,也就是用乘法原理组合。

计算公式

上面已经给出了一个递推式:

H_n=\sum_{i=0}^{n-1} H_{i-1} H_{n-i}

还有一个递推式:

H_n=\dfrac{H_{n-1}(4n-2)}{n+1}

通项公式:

H_n=\dfrac{C_{2n}^n}{n+1}=C_{2n}^n-C_{2n}^{n-1}

使用递推式可以在 O(n) 复杂度内预处理 H_i。而使用组合数通项公式再用 lucas / exLucas 等科技则可以获得更低的复杂度。

Stirling 数

第二类 Stirling 数

第二类 Stirling 数 S(n,k) 表示将 n 个不同元素分为 k 个非空子集的方案数。

第二类 Stirling 数有递推计算和通项计算两种计算方法。

递推式

考虑到第 n 个元素时,有两种方案:

  1. 将第 n 个元素单独放一个新子集,方案数为 S(n-1,k-1)
  2. 将第 n 个元素放一个已有子集,方案数为 k \times S(n-1,k)

由加法原理得:

S(n,k)=S(n-1,k-1)+k \times S(n-1,k)

通项公式

S(n,k)=\sum_{i=0}^{k} \dfrac{(-1)^{k-i} i^n}{i!(k-i)!}

证明不会。

第一类 Stirling 数

第一类 Stirling 数 s(n,k) 表示将 n 个不同元素分为 k 个非空环形排列的方案数。

第一类 Stirling 数可通过递推计算。考虑到第 n 个元素时,仍有两种方案:

  1. 将第 n 个元素单独放一个新环,方案数为 s(n-1,k-1)
  2. 将第 n 个元素放一个已有环,方案数为 (n-1) \times s(n-1,k)

由加法原理得:

s(n,k)=s(n-1,k-1)+(n-1) \times s(n-1,k)

Bell 数

Bell 数 B_n 表示大小为 n 的集合的划分方法数。例如,对于集合 \{1,2,3\},有 5 种划分方法:

{1},{2},{3}
{1,2},{3}
{1,3},{2}
{1},{2,3}
{1,2,3}

所以 B_3=5

递推式

B_n 对应集合 \{a_1,a_2,a_3,\cdots,a_n\}B_{n+1} 对应集合 \{a_1,a_2,a_3,\cdots,a_n,a_{n+1}\},只需考虑元素 a_{n+1} 的贡献。

当它和 k 个元素分到一个子集时,还剩下 n-k 个元素,则多出的方案数有 C_{n}^{n-k} B_{n-k},且 k \in [0,n]。因此:

B_{n+1}=\sum_{k=0}^{n} C_n^k B_k

仿照杨辉三角,可以构造一个贝尔三角形:

此时 B_n=a_{n,0}。代码如下:

int b[N][N];
void pre() {
  b[0][0] = 1;
  for (int i = 1; i < N; i++) {
    b[i][0] = b[i - 1][i - 1];
    for (int j = 1; j <= i; j++) b[i][j] = b[i - 1][j - 1] + b[i][j - 1];
  }
}

与第二类 Stirling 数的关系

枚举划分成 k 个非空集合,则每种情况的方案数为第二类 Stirling 数 S(n,k)。于是:

B_n=\sum_{k=0}^n S(n,k)

这表明:Bell 数 B_n 就是第 n 行第二类 Stirling 数的和。

Eulerian 数

在组合数学中,Eulerian 数 A(n,m) 是指 1n 的排列中满足恰好 m 个元素大于前一个元素的排列数目。下面我们称一个元素大于前一个元素的位置为上升位置。

Eulerian 数可通过递推计算。考虑第 n 个元素加入带来的贡献。

  1. 将它插入到一个上升位置中或排列之前,没有贡献;
  2. A(n-1,m-1) 转移时,不能满足条件 1,共有 n-(m-1)=n-m 种方案,且贡献为 1
  3. A(n-1,m) 转移时,必须满足条件 1 才有贡献,方案数为 m+1

因此:

A(n,m)=(n-m)\times A(n-1,m-1)+(m+1) \times A(n-1,m)

边界为 A(n,0)=1A(0,m)=0

代码实现:

int A(int n, int m) {
  if (m >= n || n == 0) return 0;
  if (m == 0) return 1;
  return (n - m) * A(n - 1, m - 1) + (m + 1) * A(n - 1, m); 
}