【求组合数】各种组合数卡特兰数

最新推荐文章于 2024-08-03 21:32:38 发布

HuaweiMian_

最新推荐文章于 2024-08-03 21:32:38 发布

阅读量203

点赞数

分类专栏：基础算法与基础数据结构文章标签：算法 c++

本文链接：https://blog.csdn.net/qq_73176808/article/details/132199475

版权

基础算法与基础数据结构专栏收录该内容

24 篇文章 5 订阅

订阅专栏

很多题目需要求组合数 $C^b_a$ ，而面对题目的不同要求，求组合数的方法略有不同。下面来看看几种主要的求组合数的算法。

1、根据定义求

我们拿最基本的求法来引入。组合数的公式： $C^b_a=\frac{a!}{(a-b)!b!}$ 。因此只要根据公式即可写出。规定 $C^0_n=1$ 。给出代码来求解 $C^b_a$ ：

int C(int a, int b)
{
	if (b > a) return 0;

	int res = 1;
	for (int i = 1, j = a; i <= b; i ++ , j -- )
	{
		res = res * j;
		res = res / i;
	}
	return res;
}

如果数据范围略大，需要模上一个数，这时候除法使用逆元来求 $C^a_bmodp$ ：

int C(int a, int b, int p)
{
	if (b > a) return 0;

	int res = 1;
	for (int i = 1, j = a; i <= b; i ++ , j -- )
	{
		res = (long long)res * j % p;
		//inv(i)是i的逆元
		//根据需求不同，考虑使用扩展欧几里得或快速幂
		res = (long long)res * inv(i) % p;
	}
	return res;

2、预处理结果（数的范围较小）

首先来看一道例题：
原题acwing 求组合数 I

给定 n 组询问，每组询问给定两个整数 $a$ ， $b$ ，请你输出 $C^b_amod(10^9+7)$ 的值。
输入格式
第一行包含整数 n。
接下来 n 行，每行包含一组 a 和 b。
输出格式
共 n 行，每行输出一个询问的解。
数据范围
1≤n≤100000,
1≤b≤a≤2000
输入样例：

3
3 1
5 3
2 2

输出样例：

3
10
1

处理每个组合数的时间为m，总的时间复杂度就是 $O (mn)$ ，大概是 $2×10^9$ ，明显会TLE。这时候就要用一个递推公式进行预处理了： $C^b_a=C^b_{a-1}+C^{b-1}_{a-1}$ 。这个公式的理解是：可以把所有组合分为两种，一种是不选第b个，一种是选第b个，这样就能不重不漏的包括所有组合。
有了这个公式，只要预处理每个组合数即可。这样时间复杂度就是 $O(m^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 ++ )
            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;
    cin >> n;
    
    while (n -- )
    {
        int a, b;
        scanf("%d%d", &a, &b);
        
        printf("%d\n", c[a][b]);
    }
    
    return 0;
}

3、预处理阶乘（数的范围较大）

先来看一道例题：
原题acwing 求组合数 II

给定 n 组询问，每组询问给定两个整数 $a$ ， $b$ ，请你输出 $C^b_amod(10^9+7)$ 的值。
输入格式
第一行包含整数 n。
接下来 n 行，每行包含一组 a 和 b。
输出格式
共 n 行，每行输出一个询问的解。
数据范围
1≤n≤10000,
1≤b≤a≤10⁵
输入样例：

3
3 1
5 3
2 2

输出样例：

3
10
1

可以看到题面基本和上一道题没什么差异，但是询问变为一万次，数据范围变成了十万。这样如果再用上面的做法，时间复杂度为 $O(m^2)$ ，是10¹⁰，也明显TLE了。根据公式： $C^b_a=\frac{a!}{(a-b)!b!}$ ，发现可以预处理所有阶乘，注意除法使用逆元来完成，所以也要预处理一下阶乘的逆元。这样时间复杂度就直接变成O(m)了，可以通过。给出代码：

#include <iostream>

using namespace std;

typedef long long LL;

const int N = 1e5 + 10, mod = 1e9 + 7;

int fact[N], infact[N];

LL qmi(int a, int k, int p)
{
    LL res = 1;
    while (k)
    {
        if (k & 1) res = res * a % p;
        a = (LL)a * a % p;
        k >>= 1;
    }
    
    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;
    }
    
    int n;
    cin >> 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;
}

4、卢卡斯（Lucas）定理（数的范围巨大）

先来看一道例题：
原题acwing 求组合数 III

给定 n 组询问，每组询问给定三个整数 a,b,p，其中 p 是质数，请你输出 $C^b_amodp$ 的值。
输入格式
第一行包含整数 n。
接下来 n 行，每行包含一组 a,b,p。
输出格式
共 n 行，每行输出一个询问的解。
数据范围
1≤n≤20,
1≤b≤a≤10¹⁸,
1≤p≤10⁵,
输入样例：
3
5 3 7
3 1 5
6 4 13
输出样例：
3
3
2

可以看到，这道题目的数的范围达到了10¹⁸，如果用之前的做法，那是肯定不行的。一方面，时间复杂度上不允许；另一方面，long long类型的范围大概到10¹⁹，两个数一乘就直接溢出了。这里就要介绍一种新的做法：卢卡斯定理。直接给出公式：
$C^b_a≡C^{b modp}_{amodp}·C^{⌊b/p⌋}_{⌊a/p⌋}(modp)$ ， $p$ 是质数

根据公式，我们能直接写出代码：

#include <iostream>

using namespace std;

typedef long long LL;

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

int C(int a, int b, int p)
{
    if (b > a) return 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) % p;
    }
    
    return res;
}

int lucas(LL a, LL b, int p)
{
    if (a < p && b < p) return C(a, b, p);
    return (LL)C(a % p, b % p, p) * lucas(a / p, b / p, p) % p;
}

int main()
{
    int n;
    cin >> n;
    
    while (n -- )
    {
        LL a, b;
        int p;
        cin >> a >> b >> p;
        
        cout << lucas(a, b, p) << endl;
    }
    
    return 0;
}

下面给出一个定理的证明，这个不需要掌握。

卢卡斯定理的内容为：对于非负整数 $a$ 、 $b$ ，和质数 $p$ ，有 $C^b_a≡∏^k_{i=0}C^{b_i}_{a_i}(modp)$ ，其中 $a=a_kp^k+a_{k-1}p^{k-1}+…+a_0p^0$ ， $b=b_kp^k+b_{k-1}p^{k-1}+…+b_0p^0$ 是 $a, b$ 的 $p$ 进制展开。

证明：
对于 $1+x)^p$ ，根据二项式定理，有
$(1+x)^p=C^0_px^0+C^1_px^1+…+C^p_px^p$ ；

而 $(1+x)^p=C^0_px^0+C^1_px^1+…+C^p_px^p≡C^0_px^0+C^p_px^p（modp)$ ，因为 $p$ 是质数；

也即 $1+x)^p≡1+x^p$ ， $p$ 是质数；

对于 $1+x)^a$ ，令 $a=a_kp^k+a_{k-1}p^{k-1}+…+a_0p^0$ ，
有 $1+x)^a=((1+x)^{p^0})^{a_0}((1+x)^{p^1})^{a_1}…((1+x)^{p^k})^{a_k}$ ；
则 $1+x)^a≡(1+x^{p^0})^{a_0}(1+x^{p^1})^{a_1}…(1+x^{p^k})^{a_k}(modp)$

对比 $x^b$ 的系数，对于左侧，系数为 $C^b_a$ ，对于右侧，令 $b=b_kp^k+b_{k-1}p^{k-1}+…+b_0p^0$ ，那么系数为 $C^{b_0}_{a_0}C^{b_1}_{a_1}…C^{b_k}_{a_k}$ ；

因此有 $C^b_a≡C^{b_0}_{a_0}C^{b_1}_{a_1}…C^{b_k}_{a_k}(modp)$ ，得证。

另外，如果 $b > a$ ，那么右侧是不存在这样的系数的。所以当 $b > a$ 时， $C^b_a=0$ 。

对于我们用到的公式，只要结合卢卡斯定理的证明以及 $a=⌊\frac{a}{p}⌋p+amodp$ ，即可得出结论。这里不再赘述。

5、不模的组合数

先来看例题：
原题acwing 求组合数 IV

输入 a,b，求 $C^b_a$ 的值。
注意结果可能很大，需要使用高精度计算。
输入格式
共一行，包含两个整数 a 和 b。
输出格式
共一行，输出 $C^b_a$ 的值。
数据范围
1≤b≤a≤5000
输入样例：
5 3
输出样例：
10

这个问题就要从定义出发了，不模的话，结果的范围很容易就溢出了。这时候就要使用高精度来进行计算了。这里使用组合数公式： $C^b_a=\frac{a!}{(a-b)!b!}$ 。但是这里还涉及到除法，所以我们选择把整个数分解质因数，然后把质因数乘在一起即可。

首先要关注怎么处理除法。对于除以一个数，只需要把它质因数的幂次减掉即可。
其次关注怎么求 $n!$ 的质因数 $p$ 的幂次。给出公式：
$n!$ 中质因数 $p$ 的幂次 $=⌊\frac{a}{p}⌋+⌊\frac{a}{p^2}⌋+⌊\frac{a}{p^3}⌋+…$
证明也比较容易，只需要依次加上 $p$ 的倍数的个数、 $p^2$ 的倍数的个数……即可。

给出代码：

#include <iostream>
#include <vector>

using namespace std;

const int N = 5010;

int primes[N], cnt;
int sum[N];
bool st[N];

//筛质数
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;
        }
    }
}

//求n的阶乘中质因数p的幂次
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() || 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;
    
    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 -- ) printf("%d", res[i]);
    puts("");
    
    return 0;
}

6、组合数应用：卡特兰数

先由一道例题引入
原题acwing 满足条件的01序列

给定 n 个 0 和 n 个 1，它们将按照某种顺序排成长度为 2n 的序列，求它们能排列成的所有序列中，能够满足任意前缀序列中 0 的个数都不少于 1 的个数的序列有多少个。
输出的答案对 10⁹+7 取模。
输入格式
共一行，包含整数 n。
输出格式
共一行，包含一个整数，表示答案。
数据范围
1≤n≤10⁵
输入样例：
3
输出样例：
5

首先把这道题目抽象为：给定坐标轴，要走到 $(n, n)$ 点，序列中0代表向左走一格，1代表向上走一格。那么符合要求的路线就是对于走到的每个点，向左走的距离一定大于等于向上走的距离。抽象到坐标轴上，即路线一定要在红线之下，也即路线一定不能经过红线上的任何一点。
在这里插入图片描述
这里假设走到 $(6, 6)$ 。容易得出，走到 $(6, 6)$ 的所有走法数量为 $C^6_{12}$ （十二步里选六步向左走），现在要做的就是用所有走法减去不符合要求的走法。随便举一个不符合要求的例子：

处理方法为：找到第一个与红线的交点，并将此点后的路线关于红线做对称。得到：
在这里插入图片描述
不难得到：所有不符合条件的到 $(6, 6)$ 的路线（蓝色）都可以通过做对称化为一条到 $(5, 7)$ 的路线（绿色），并且它们是一一对应的。因此，不符合条件的路线数就是 $C^5_{12}$ 。那么结果即为 $C^6_{12}-C^5_{12}$ 。

进行推广，假设要走到 $(n, n)$ ，那么合法路线数就是 $C^n_{2n}-C^{n-1}_{2n}$ 。
对于满足此性质的所有序列的个数，就叫做卡特兰数。它是组合数学中一个重要的结论，很多算法题里都要用到卡特兰数。

现在对公式进行变形：
$C^n_{2n}-C^{n-1}_{2n}$
$=\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}\frac{2n!}{n!n!}$
$=\frac{1}{n+1}C^n_{2n}$
得到最终的公式： $C^n_{2n}-C^{n-1}_{2n}=\frac{C^n_{2n}}{n+1}$ 。

有了这个，即可给出代码：

#include <iostream>

using namespace std;

typedef long long LL;

const int mod = 1e9 + 7;

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

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