算法学习笔记（40）——组合计数

组合计数

• 组合计数
  • 求组合数 I —— 预处理组合数
  • 求组合数 II —— 预处理阶乘
  • 求组合数 III —— 卢卡斯定理（Lucas）
  • 求组合数 IV —— 高精度
  • 满足条件的01序列 —— 卡特兰数

求组合数 I —— 预处理组合数

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

数据范围：
$1 \le n \le 10000$
$1 \le b \le a \le 2000$

如果直接暴力做的话，每次计算一个组合数，最坏需要2000次运算（因为a与b的范围是2000），而最多有10000次询问，所以时间复杂度达到2千万（ $2000 * 10000 = 2 * 10^7$ ），会超时，所以不能这么算。
我们注意到 $a$ 和 $b$ 的范围是2000，所以不同的 $a$ 和 $b$ 的方案数有 $2000^2$ 组，总共只有四百万对（ $2000^2 = 4 * 10^6$ ），所以我们可以先预处理出所有方案的情况，利用以下递推式：

$$C_a^b = C_{a-1}^b + C_{a-1}^{b-1}$$

$C_a^b$ 表示从 $a$ 个苹果中选 $b$ 个的方案数，假设 $a$ 个苹果中有一个红色苹果，我们可以将方案分为两大类：

• 第一类的方案中选择的 $b$ 个苹果包含红色苹果，则这种情况下的方案数就是除了红色苹果之外，从剩下的苹果中再选择 $b-1$ 个苹果： $C_{a-1}^{b-1}$
• 第二类的方案中选择的 $b$ 个苹果不包含红色苹果，则这种情况下的方案数就是在除去红色苹果的剩余部分选择 $b$ 个苹果：$C_{a-1}^b$

题目链接：AcWing 885. 求组合数 I

时间复杂度：$O(N^2)$

#include <iostream>

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 ++ )
            // 处理边界情况，从i个苹果里选0个苹果的方案数是1
            if (!j) c[i][j] = 1;
            else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;
}

int main()
{
    init();
    
    int n;
    scanf("%d", &n);
    while (n -- ) {
        int a, b;
        scanf("%d%d", &a, &b);
        printf("%d\n", c[a][b]);
    }
    
    return 0;
}

求组合数 II —— 预处理阶乘

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

数据范围：
$1 \le n \le 10000$
$1 \le b \le a \le 10^5$

此时如果用上一道题的预处理方法，则由于 $a$ 和 $b$ 的范围是十万，不同的 $a$ 和 $b$ 的方案数有一百亿组（ $10^{10}$ ），超时QAQ。所以我们需要换一种预处理方式。我们知道组合数公式为：

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

所以可以将公式中的阶乘预处理出来，$fact[i] = i! \bmod (10^9+7)$
又因为

$$\frac{a}{b} \bmod p \neq \frac{a \bmod p}{b \bmod p}$$

可以通过求逆元 $x$（ $\frac{a}{b} \bmod p = a x \bmod p$ ），将除法转换为乘法，再预处理出 $infact[i] = (i!)^{-1} \bmod (10^9+7)$，则

$$C_n^m = fact[n] \times infact[n-m] \times infact[m]$$

本题中模数是质数，符合费马小定理的条件，可以使用快速幂求逆元。
快速幂的时间复杂度为$O(\log k)$，$k$ 是幂次数（我们将幂次分为 $2^0, 2^1,2^2 \cdots$）

题目链接：AcWing 886. 求组合数 II

时间复杂度： $O(N \log N)$

#include <iostream>

using namespace std;

typedef long long LL;

const int N = 100010;
const int mod = 1e9 + 7;

int n;
LL fact[N], infact[N];

// 快速幂，返回a^k % p
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()
{
    // 预处理阶乘数组
    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;
    }
    
    scanf("%d", &n);
    while (n -- ) {
        int a, b;
        scanf("%d%d", &a, &b);
        // 两个10^9相乘不会溢出long long，三个相乘有可能溢出，所以中间取模一次
        printf("%d\n", fact[a] * infact[b] % mod * infact[a - b] % mod);
    }
    
    return 0;
}

求组合数 III —— 卢卡斯定理（Lucas）

给定 n 组询问，每组询问给定三个整数 $a,b,p$，其中 $p$ 是质数，请你输出 $C_b^a\bmod p$ 的值。

数据范围：
$1 \le n \le 20$
$1 \le b \le a \le 10^{18}$
$1 \le p \le 10^5$

Lucas 定理
若 $p$ 是质数，则对于任意整数 $1 \le m \le n$ ，有:

$$C_n^m \equiv C_{n \bmod p}^{m \bmod p} * C_{n / p}^{m / p} \pmod p$$

也就是把 $n$ 和 $m$ 表示成 $p$ 进制数，对 $p$ 进制下的每一位分别计算组合数，最后再乘起来。证明一般需要用到生成函数知识，此处仅简略证明。

证明：
把 $n$ 和 $m$ 表示成 $p$ 进制数

$$\begin{align*} & n = n_kp^k+ n_{k-1}p^{k-1} + \dots + n_0p^0 \\ & m = m_kp^k+ m_{k-1}p^{k-1} + \dots + m_0p^0 \end{align*}$$

由二项式定理 $(a+b)^n = \sum\limits_{k=0}^{n} C_n^ka^kb^{n-k}$ 可以得到：

$$\begin{align*} (1+x)^p & = C_p^0 \cdot x^0 \cdot 1^p + C_p^1 \cdot x^1 \cdot 1^{p-1} + \dots + C_p^p \cdot x^p \cdot 1^0 \\ & = C_p^0 \cdot 1 + C_p^1 \cdot x + \dots + C_p^p x^p \end{align*}$$

上式第1项到第n-1项分子都有 $p$ ，对 $p$ 取模都得零，于是：

$$(1+x)^p \equiv 1+x^p \pmod p$$

将 $n$ 分解为 $p$ 进制数，我们得到：

$$\begin{align*} (1+x)^n & = (1+x)^{n_0p^0 + n_1p^1 + n_2p^2 + \dots + n_kp^k} \\ & = (1+x)^{n_0} \cdot ((1+x)^p)^{n_1} + ((1+x)^{p^2})^{n_2} + \dots + ((1+x)^{p^k})^{n_k} \\ & \equiv (1+x)^{n_0} \cdot (1+x^{p})^{n_1} \cdot (1+x^{p^2})^{n_2} + \dots + (1+x^{p_k})^{n_k} \end{align*}$$

二项式定理展开，上式左边 $x^m$ 的系数是 $C_n^m$ ，右边的系数是 $C_{n_k}^{m_k} \cdot C_{n_{k-1}}^{m_{k-1}} \dots C_{n_0}^{m_0}$，则

$$C_n^m \equiv C_{n_k}^{m_k} \cdot C_{n_{k-1}}^{m_{k-1}} \dots C_{n_0}^{m_0} \pmod p$$

卢卡斯定理就是上式的递归形式。

证毕

题目链接：AcWing 887. 求组合数 III

将 $n$ 和 $m$ 表示成 $p$ 进制数的时间复杂度是 $O(\log_p 10^{18})$ ，利用快速幂从定义出发（阶乘除法）计算一个 $10^5$ 大小的组合数的时间复杂度是 $O(p \log p)$，而当 $p$ 取极限值时，$\log_p 10^{18}$ 很小，所以总的时间复杂度约为：$O(p \log p)$

#include <iostream>

using namespace std;

typedef long long LL;

int p;

// 快速幂，返回a^k % p，p在本题是全局变量，所以不需要作为参数引入
LL qmi(int a, int k)
{
    LL res = 1;
    while (k) {
        if (k & 1) res = (LL)res * a % p;
        k >>= 1;
        a = (LL)a * a % p;
    }
    return res;
}

// 返回组合数C_a^b = a! / b! * (a-b)!
int C(int a, int b)
{
    int res = 1;
    for (int i = 1, j = a; i <= b; i ++, j -- ) {
        // 组合数中a!与(a-b)!约分，则分子只需累乘a到a-b+1
        res = (LL)res * j % p;
        // 分母累积b! 等价于乘以对应的逆元
        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;
}

求组合数 IV —— 高精度

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

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

数据范围：
$1 \le b \le a \le 5000$

如果从定义出发，我们需要使用高精度乘法和高精度除法，代码实现比较复杂，所以我们首先对组合数分解质因数，从而只需用高精度乘法即可计算结果。

$$C_a^b = \frac{a!}{b!(a-b)!} = p_1^{\alpha_1} \cdot p_2^{\alpha_2} \cdot p_3^{\alpha_3} + \dots$$

而 $a!$ 内包含 $p$ 的次数为：

$$\lfloor \frac{a}{p} \rfloor + \lfloor \frac{a}{p^2} \rfloor + \lfloor \frac{a}{p^3} \rfloor + \dots$$

$\lfloor \frac{a}{p} \rfloor$ 表示 $1$ ～ $a$ 中 $p$ 的倍数有多少个
$\lfloor \frac{a}{p^2} \rfloor$ 表示 $1$ ～ $a$ 中 $p^2$ 的倍数有多少个
以此类推...

举例说明，$1,2,3,4,5,6,7,8$ 中，在 $2,4,6,8$ 中均出现了4次 $2$ ，在 $4,8$ 中出现了2次 $2^2$ ，在 $8$ 中出现了1次 $2^3$，于是 $8!$ 中共包含7个 $2$。
该步骤时间复杂度是 $O(\log p)$

算法思路：

1. 筛质数，将1～5000内的质数筛出来
2. 求每个质数的次数
3. 用高精度乘法把所有的质因子乘到一起

#include <iostream>
#include <vector>

using namespace std;

const int N = 5010;

int primes[N], cnt; // 存储筛选出的质数与对应的数组下标
bool st[N];         // 标记是质数还是合数
int sum[N];         // 存储阶乘中包含的每个质数的次数

// 线性筛质数
void get_primes(int n)
{
    // 质数从2开始
    for (int i = 2; i <= n; i ++ ) {
        if (!st[i]) primes[cnt ++] = i;
        for (int j = 0; primes[j] <= n / i; j ++ ) {
            st[primes[j] * i] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

// 返回a!里面包含的p的次数
int get(int a, int p)
{
    int res = 0;
    while (a) {
        res += a / p;   // a/p + a/p^2 + a/p^3 + ...
        a /= p;         // 当a < p^k时, a = a / p^k = 0, 跳出循环
    }
    return res;
}

// 高精度乘法模板
vector<int> mul(vector<int> a, int &b)
{
   vector<int> c;
   int t = 0;
   for (int i = 0; i < a.size() || t; i ++ ) {
       if (i < a.size()) t += a[i] * b;
       c.push_back(t % 10);
       t /= 10;
   }
   while (c.size() > 1 && c.back() == 0) c.pop_back();
   return c;
}

int main()
{
    int a, b;
    cin >> a >> b;
    
    // 筛出1～N范围内的质数
    get_primes(a);
    
    // 求每个质数的次数
    for (int i = 0; i < cnt; i ++ ) {
        int p = primes[i];
        sum[i] = get(a, p) - get(b, p) - get(a - 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;
}

满足条件的01序列 —— 卡特兰数

题目链接:AcWing 889. 满足条件的01序列

01排列的方案数与从(0,0)走到(6,6)的路径数一致。任意前缀中0的个数都要小于1的个数，路径上的点的横纵坐标需满足不经过下图红色线上的点。

总共有12步，其中6步向上走，所以从(0,0)走到(6,6)的所有路径数是 $C_{12}^6$ 。
任何一条从(0,0)走到(6,6)并且经过红边的路径数量（从第一个与红边相交的点做轴对称，可以对应一条从(0,0)走到(5,7)的路径），从而等于从(0,0)走到(5,7)的路径数量 $C_{12}^5$ 。
所以满足条件的路径数为：$C_{12}^6 - C_{12}^5$
通用形式为：

$$\begin{align*} C_{2n}^{n} - C_{2n}^{n-1} &= \frac{(2n)!}{n!n!} - \frac{(2n)!}{(n-1)!(n+1)!} \\ &= \frac{(2n)!(n+1) - (2n)!n}{(n+1)!n!} \\ &= \frac{(2n)!}{(n+1)!n!} \\ &= \frac{1}{n+1} \cdot \frac{(2n)!}{n!(2n-n)!} \\ &= \frac{C_{2n}^{n}}{n+1} \end{align*}$$

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

$$Cat_n = \frac{C_{2n}^{n}}{n+1}$$

#include <iostream>

using namespace std;

typedef long long LL;

const int mod = 1e9 + 7;

// 模数为质数，所以可以根据费马小定理，用快速幂求逆元
int qmi(int a, int k, int p)
{
    int res = 1;
    while (k) {
        if (k & 1) res = (LL)res * a % mod;
        a = (LL)a * a % mod;
        k >>= 1;
    }
    return res;
}

int main()
{
    int n;
    cin >> n;
    
    int a = 2 * n, b = n;
    int res = 1;
    
    // (2n)! / n!
    for (int i = a; i > a - b; i -- ) res = (LL)res * i % mod;
    // 1 / n!
    for (int i = 1; i <= b; i ++ ) res = (LL)res * qmi(i, mod - 2, mod) % mod;
    // 1 / n + 1
    res = (LL)res * qmi(n + 1, mod - 2, mod) % mod;
    
    cout << res << endl;
    
    return 0;
}