组合数

组合数

1. 组合数原理

组合数

• 这里的组合数指的是排列组合中的组合数。
• 组合中的问题是：给定a个苹果，从中选择出b个苹果，请问有多少种选法？结果就是组合数$C_a^b$。
• 关于 $C_a^b$ 的求解公式如下：

$$C_a^b=\frac{a!}{b!\times (a-b)!}=\frac{a\times (a-1)\times ...\times (a-b+1)}{1\times 2\times ...\times b}$$

• 另外一个关于 $C_a^b$ 的递推公式：

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

从实际含义考虑这个式子即可：从a个苹果中选择出b个苹果 = （首先拿出一个苹果）从a-1个苹果中选择出b-1个苹果（包含最开始拿出的一个苹果） + 从a-1个苹果中选择出b个苹果（不包含最开始拿出的一个苹果）

• 下面的四个组合数的题目主要是数据范围不一样，如下：

/**
 * 求n组数据的组合数C(a, b)
 * AcWing 885. 求组合数 I
 * 		1 <= n <= 10000, 1 <= b <= a <= 2000, p = 10^9+7
 * AcWing 886. 求组合数 II
 * 		1 <= n <= 10000, 1 <= b <= a <= 10^5, p = 10^9+7
 * AcWing 887. 求组合数 III
 * 		1 <= n <= 20, 1 <= b <= a <= 10^18, 1 <= p <= 10^5且p是质数
 * 
 * AcWing 888. 求组合数 IV
 *   	求出具体的组合数的值  1 <= b <= a <= 5000
 */

• 关于不同数据范围组合数的求解具体方法请参考下面的题目分析。

---

卡特兰数

• 卡特兰数：首先这个一个数，很多问题的结果都是卡塔兰数，比如2016年全国三卷数学选择题压轴题让求解的就是卡特兰数，问题如下（AcWing 889. 满足条件的01序列也是这个题）：

• 首先是结论：卡特兰数为:

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

• 因此，对于上面的题目，结果就是

$$\frac{C_{2m}^m}{m+1}=\frac{C_8^4}{4+1}=\frac{70}{5}=14$$

因此选择C。

• 下面看一下卡特兰数的公式是如何推导出来的。首先我们需要将上述问题转换成一个等价的问题：在一个二维平面内，从(0, 0)出发到达(n, n)，每次可以向上或者向右走一格，0代表向右走一个，1代表向上走一格，则每条路径都会代表一个01序列，则满足任意前缀中0的个数不少于1个数序列对应的路径则右下侧，如下图：

符合要求的路径必须严格在上图中红色线的下面（不可以碰到图中的红线，可以碰到绿线）。则我们考虑任意一条不合法路径，例如下图：

所有路径的条数为 $C_{2n}^n$，其中不合法的路径有 $C_{2n}^{n-1}$ 条，因此合法路径有：
$$\begin{gathered} C_{2n}^n-C_{2n}^{n-1}=\frac{(2n)!}{n!\times n!}-\frac{(2n)!}{(n+1)!\times (n-1)!} \\ =\frac{(2n)!\times (n+1)}{n!\times (n+1)!}-\frac{(2n)!\times n}{n!\times (n+1)!} \\ =\frac{(2n)!\times (n+1)-(2n)!\times n}{n!\times (n+1)!} \\ =\frac{(2n)!}{n!\times (n+1)!}=\frac{1}{n+1}\times \frac{(2n)!}{n!\times n!} \\ =\frac{C_{2n}^n}{n+1} \end{gathered}$$
推导完毕。

• 除了上述两种问题，如下问题对应的答案也是卡特兰数：

  （1）n个结点的二叉树数量h(n) ；其实有递推公式，即：
  $$h(n)=\sum _{i=1}^{n}h(i-1)\times h(n-i)h(0)=1$$
  （2）矩阵链乘：$P=A_1\times A_2\times ...\times A_n$，有多少种不同的计算次序？（相当于加括号，问合法括号序列有多少个）

  （3）一个栈(无穷大)的进栈序列为1，2，3，…，n，有多少个不同的出栈序列?

  （4）有2n个人排成一行进入剧场。入场费5元。其中只有n个人有一张5元钞票，另外n人只有10元钞票，剧院无其它钞票，问有多少中方法使得只要有10元的人买票，售票处就有5元的钞票找零？(将持5元者到达视作将5元入栈，持10元者到达视作使栈中某5元出栈)

• 补充内容：

2. AcWing上的组合数题目

AcWing 885. 求组合数 I

问题描述

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

  

分析

• 这一题可以使用组合数的递推公式求解，因为$1\le b\le a\le 2000$，计算量最大为$2000^2=4\times 10^6$，大约四百万的复杂度，完全可以接受。

• 递推公式：$C_a^b=C_{a-1}^{b-1}+C_{a-1}^b$。

代码

• C++

#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++)
            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;
}

AcWing 886. 求组合数 II

问题描述

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

  

分析

• 这一题如果仍然使用递推的方式求解的话，因为$1\le b\le a\le 10^5$，计算量最大为$10^{10}$，大约十亿的复杂度，完全不能被接受。

• 这一题需要使用到的公式是：

$$C_a^b=\frac{a!}{b!\times (a-b)!}$$

• 存在的问题是阶乘太大，计算机无法存储下来，另外我们最终的结果需要对$10^9+7$取模，因此考虑值存取取模后的结果，但是遗憾的是：

$$\frac{a}{b}(modp)\ne \frac{a(modp)}{b(modp)}$$

• 那该怎么办呢？这里就要使用到了逆元，关于逆元可以参考：网址。

• $\frac{a}{b}\equiv a\times x(modm)$，我们求出分母的逆元即可，因为p是质数，所以任何小于p的数都与p互质，因此逆元存在，另外由于p是质数，可以使用快速幂求解逆元。

• 另外根据逆元定义很容易验证b、c的逆元相乘的结果等于$b\times c$的逆元。

• 实际代码中我们可以使用数组fact记录各个数据阶乘的结果，数组infact记录阶乘对应的逆元，最终

$$C_a^b=((fact[a]\times infact[b])\%p)\times infact[a-b])\%p$$

• 计算量：预处理各个数的阶乘以及阶乘的逆元，求逆元可以使用快速幂，因此时间复杂度为$O(n\times log(n))$，计算量为$10^5\times log(10^5)\approx 1.6\times 10^6$。可以接受。

代码

• C++

#include <iostream>

using namespace std;

typedef long long LL;

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

int fact[N];  // fact[i]存储 (i!)%mod 的结果
int infact[N];  // infact[i]存储 (i!模mod的逆元)%mod 的结果

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

int main() {
    
    fact[0] = infact[0] = 1;  // 0!=1, 1的逆元是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;
    }
    
    int n;
    scanf("%d", &n);
    while (n--) {
        int a, b;
        scanf("%d%d", &a, &b);
        printf("%d\n", (LL)fact[a] * infact[b] % mod * infact[a - b] % mod);
    }
    
    return 0;
}

AcWing 887. 求组合数 III

问题描述

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

  

分析

• 上两题中的方法这里都不能使用，因为b、a范围太大了。本题需要使用到卢卡斯定理(lucas)。

• 卢卡斯定理（要求式中a、b为非负整数，p为质数）：

$$C_a^b=C_{a(modp)}^{b(modp)}\times C_{a/p}^{b/p}(modp)$$

递归求解即可。

• 关于$C_{a(modp)}^{b(modp)}$的求解，直接使用定义求解即可：

$$C_a^b=\frac{a\times (a-1)\times ...\times (a-b+1)}{1\times 2\times ...\times b}$$

• 设$N=10^{18}$，求解$C_{a(modp)}^{b(modp)}$的复杂度为$O(p\times log(p))$的（因为要求解 $2^{p-2}\%p$ 的值，然后预处理所有的阶乘以及逆元的阶乘），后面的 $C_{a/p}^{b/p}$ 可以递归求解，可以看成p进制的数据，因此需要递归$lo{g}_p(N)$次，因此总时间复杂度为：$lo{g}_p(N)\times p\times log(p)$。一般来说p值比较大，因此$lo{g}_p(N)$比较小，可以忽略，时间复杂度是$O(p\times log(p))$的，和上一题的时间复杂度差不多。

• 关于卢卡斯定理的证明，不要求掌握。

代码

• C++

#include <iostream>

using namespace std;

typedef long long LL;

int p;

// 计算 a^b % p
int qmi(int a, int b) {
    
    int res = 1 % p;
    while (b) {
        if (b & 1) res = (LL)res * a % p;
        a = (LL)a * a % p;
        b >>= 1;
    }
    return res;
}

// 根据定义计算C(a, b) % p
int C(int a, int b) {
    
    if (b > a) return 0;  // 不合法情况，直接返回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

问题描述

• 问题链接：AcWing 888. 求组合数 IV

  

分析

• 这一题根据定义计算即可：

$$C_a^b=\frac{a!}{b!\times (a-b)!}$$

• 我们需要实现一个高精度乘法和一个高精度除法。这样的写法运行效率比较低，另外也比较难写。
• 我们可以考虑对$C_a^b$进行质因数分解，假设：

$$C_a^b=p_1^{{\alpha }_1}\times p_2^{{\alpha }_2}\times ...\times p_k^{{\alpha }_k}$$

这样我们只需要写高精度乘法即可。

• 如何计算每个质因子出现的次数呢？我们发现$C_a^b$是一个阶乘除以另外两个阶乘，因此，我们可以分别求出这三个阶乘的质因数分解，然后将对应质因数的幂次相减即可。此时问题就转换成了如何求解阶乘的质因数分解，对应题目：AcWing 197. 阶乘分解。

• 阶乘的质因数分解：$n!$所有的质因子一定是小于等于n的，否则如果存在大于n的质因子，说明该质因子一定是某几个数相乘得到的，违反质数的定义。

  （1）求出1~n中所有的质数，可以使用线性法筛质数，可以参考：网址；

  （2）枚举某个质数p，则其出现的次数为：
  $$\lfloor \frac{n}{p}\rfloor +\lfloor \frac{n}{p^2}\rfloor +......$$
  该方法求阶乘$n!$的质因数分解的时间复杂度大约为$O(n)$，因为1~n中大约有$\frac{n}{log(n)}$个质数，求每个质数出现次数计算大约为$log(n)$的，相乘得到阶乘质因数分解时间复杂度大约为$O(n)$的。

• 本题的解题步骤是：

  （1）使用线性筛法筛出1~a之间的所有质数，时间复杂度$O(a)$；

  （2）根据 $C_a^b$ 定义计算其质因数分解结果，时间复杂度$O(a)$；

  （3）高精度乘法计算最终结果。

代码

• C++

#include <iostream>
#include <vector>

using namespace std;

const int N = 5010;

int primes[N], cnt;
bool st[N];  // 该数据是否被筛掉
int sum[N];  // sum[i]代表primes[i]在C(a, b)的质因数分解中出现的次数

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[primes[j] * i] = 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;
    
    // (1) 使用线性筛法筛出1~a之间的所有质数
    get_primes(a);
    // (2) 根据 C(a, b) 定义计算其质因数分解结果
    for (int i = 0; i < cnt; i++) {
        int p = primes[i];
        sum[i] = get(a, p) - get(b, p) - get(a - b, p);
    }
    // (3) 高精度乘法计算最终结果
    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--) printf("%d", res[i]);
    puts("");
    
    return 0;
}

AcWing 889. 满足条件的01序列

问题描述

• 问题链接：AcWing 889. 满足条件的01序列

  

分析

• 本题让求卡特兰数，直接带入公式即可：

$$\frac{C_{2n}^n}{n+1}=\frac{1}{n+1}\times \frac{a\times (a-1)\times ...\times (a-b+1)}{1\times 2\times ...\times b}a=2n,b=n$$

• 因为mod是质数，因此任何数与mod都互质，可以使用快速幂求逆元，否则需要使用扩展欧几里得算法求逆元。

代码

• C++

#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 % p;
    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 = 2 * n, 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;
}