数学知识3.1-组合数

一、简述

记录一下组合数的求解问题。

二、组合数

1.[AcWing885.求组合数I]-递推

题目描述

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

输入格式

第一行包含整数 $n$。

接下来 $n$ 行，每行包含一组 $a$ 和 $b$。

输出格式

共 $n$ 行，每行输出一个询问的解。

数据范围

$1\le n\le 10000,$

$1\le b\le a\le 2000$

输入样例

3
3 1
5 3
2 2

输出样例

3
10
1

解题思路

递推。数据范围不大，根据公式 $C_a^b=C_{a-1}^b+C_{a-1}^{b-1}$，可在 $O(n^2)$ 的时间复杂度内预处理出来，然后对于每次询问 $O(1)$ 的时间输出。

关于 $C_a^b=C_{a-1}^b+C_{a-1}^{b-1}$ 的理解：从 $a$ 个苹果中选择 $b$ 个，假设其中一个苹果是 $s$，那么我们选择 $b$ 个苹果无非是包括 $s$ 和不包括 $s$，如果包括 $s$，那么 $s$ 一定被选，剩余 $a-1$ 个里面选择 $b-1$ 个，即 $C_{a-1}^{b-1}$；否则 $s$ 不被选择，那就是从 $a-1$ 个苹果里面选择 $b$ 个，即 $C_{a-1}^b$。

C++代码

#include <bits/stdc++.h>
using namespace std;
const int N = 2010, MOD = 1e9 + 7;

int n;
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();
    cin >> n;
    while (n --)
    {
        int a, b;
        cin >> a >> b;
        cout << c[a][b] << "\n";
    }
    return 0;
}

2.[AcWing886.求组合数II]-逆元

题目描述

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

输入格式

第一行包含整数 $n$。

接下来 $n$ 行，每行包含一组 $a$ 和 $b$。

输出格式

共 $n$ 行，每行输出一个询问的解。

数据范围

$1\le n\le 10000,$
$1\le b\le a\le {10}^5$

输入样例

3
3 1
5 3
2 2

输出样例

3
10
1

解题思路

数据范围不允许递推了。

首先 $C_a^b=\frac{a!}{b!(a-b)!}$，那么我们可以预处理出来阶乘，然后按照公式计算即可，但是 $\frac{a!}{b!(a-b)!}mod(1e^9+7)\ne \frac{a!mod(1e9+7)}{b!(a-b)!mod(1e^9+7)}$，但是 $\frac{a!}{b!(a-b)!}mod(1e^9+7)=(a!mod(1e^9+7))\ast (\frac{1}{b!(a-b)!}mod(1e^9+7))mod(1e^9+7)$，而 $1e^9+7$ 是质数，则根据费马定理 ($n$ 为质数，且 $a$ 与 $n$ 互质，那么 $a^{n-1}modn=1$，$a^{n-2}$ 为 $a$ 的逆元)，据此可计算每个数的逆元的阶乘。时间复杂度为 $O(NlogMOD)$。

如果 $MOD$ 不是质数，我们就不能根据费马小定理使用快速幂来计算逆元，而应使用拓展欧几里得算法来求解，具体思路：$a$ 与 $m$ 互质，$ax\equiv 1(modm)$，而 $(km+1)\equiv 1(modm)$，则可取得整数 $k$，使得 $bx=km+1$，即 $bx-km=1$，可使用拓展欧几里得算法求解。

C++代码

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 100010, MOD = 1e9 + 7;

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

LL 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;
}

void init()
{
    fact[0] = infact[0] = 1;
    for (int i = 1; i < N; i ++)
    {
        fact[i] = fact[i - 1] * i % MOD;
        infact[i] = infact[i - 1] * qmi(i, MOD - 2, MOD) % MOD;
    }
}

int main()
{
    init();
    scanf("%d", &n);
    while (n --)
    {
        int a, b;
        scanf("%d%d", &a, &b);
        LL res = fact[a] * infact[b] % MOD * infact[a - b] % MOD;
        printf("%lld\n", res);
    }
    return 0;
}

3.[AcWing887.求组合数III]-Lucas定理

题目描述

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

输入格式

第一行包含整数 $n$。

接下来 $n$ 行，每行包含一组 $a,b,p$。

输出格式

共 $n$ 行，每行输出一个询问的解。

数据范围

$1\le n\le 20,$

$1\le b\le a\le {10}^{18},$

$1\le p\le {10}^5$

输入样例

3
5 3 7
3 1 5
6 4 13

输出样例

3
3
2

解题思路

首先介绍一下 $Lucas$ 定理。

$Lucas$ 定理：$C_a^b\equiv C_{amodp}^{bmodp}\ast C_{a/p}^{b/p}(modp)$。

证明如下：
首先，$C_p^i\equiv 0(modp)(0<i<p)$，$p$ 为质数，这个从组合数的计算公式出发即可知。

然后 $(1+x{)}^p\equiv (1+x^p)modp$。这个将左边展开，利用上面的结论即可。

那么对于 $a,b$，可得 $a=np+q,b=tp+r(q<p,r<p)$。那么 $(1+x{)}^a\equiv (1+x{)}^{np+q}\equiv ((1+x{)}^p{)}^n(1+x{)}^q\equiv (1+x^p{)}^n(1+x{)}^q\equiv (1+n{x}^p+...+C_n^{n-i}{x}^{ip}+...+x^{np})(1+qx+...+C_q^{q-i}{x}^i+...+x^q)(modp)$。
这里由于 $x$ 是未知数，那么左右两侧对应 $x^i$ 的系数相同，则左侧 $x^b$ 的系数为 $C_a^{a-b}=C_a^b=C_{np+q}^{tp+r}$，对于右侧 $x^b=x^{tp}{x}^r$，即 $C_n^{n-t}{x}^{tp}\ast C_q^{q-r}{x}^r$，即得 $C_n^{n-t}\ast C_q^{q-r}=C_n^t\ast C_q^r$。那么最终得到 $C_a^b\equiv C_n^t\ast C_q^r(modp)$，即 $C_a^b\equiv C_{amodp}^{bmodp}\ast C_{a/p}^{b/p}(modp)$。

C++代码

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;

int n;

int qmi(int a, int b, int p)
{
    int res = 1;
    while (b)
    {
        if (b & 1) res = (LL) res * a % p; 
        a = (LL) a * a % p;
        b >>= 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()
{
    scanf("%d", &n);
    while (n --)
    {
        LL a, b;
        int p;
        scanf("%lld%lld%d", &a, &b, &p);
        printf("%d\n", Lucas(a, b, p));
    }
    return 0;
}

4.[AcWing888.求组合数IV]-高精度

题目描述

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

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

输入格式

共一行，包含两个整数 $a$ 和 $b$。

输出格式

共一行，输出 $C_a^b$ 的值。

数据范围

$1\le b\le a\le 5000$

输入样例

5 3

输出样例

10

解题思路

首先引入一个知识点：阶乘的分解。

我们知道阶乘是很大的一个数，那么我们如何对其进行质因数分解呢？我们考虑用质数 $p$ 来分解 $n!$，那么 $1$~$n$ 中，最多有 $\lfloor \frac{n}{p}\rfloor$ 个数是 $p$ 的倍数，那么在阶乘中也就可以分解出 $\lfloor \frac{n}{p}\rfloor$，不过 $1$~$n$ 中有些数可能是 $p^2$ 的倍数，有 $\lfloor \frac{n}{p^2}\rfloor$ 个，以此类推，我们将这些数加起来即是 $n!$ 中 $p$ 的指数。示例代码如下：

int get(int n, int p)
{
    int res = 0;
    while (n)
    {
        res += n / p;
        n /= p;
    }
    return res;
}

那么根据 $C_a^b=\frac{a!}{b!(a-b)!}$，那么关于 $p$ 的指数就是 $get(a,p)-get(b,p)-get(a-b,p)$，然后质数我们利用线性筛得到，然后进行大数乘小数的高精度乘法。

C++代码

#include <bits/stdc++.h>
using namespace std;
const int N = 5010;

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

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> res;
    int t = 0;
    for (int i = 0; i < a.size(); i ++)
    {
        t += a[i] * b;
        res.push_back(t % 10);
        t /= 10;
    }
    while (t)
    {
        res.push_back(t % 10);
        t /= 10;
    }
    return res;
}

void get_primes(int n)
{
    for (int i = 2; i <= n; i ++)
    {
        if (!st[i]) primes[cnt ++] = i;
        for (int j = 0; primes[j] * i <= n; j ++)
        {
            st[primes[j] * i] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

int main()
{
    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> ans;
    ans.push_back(1);
    for (int i = 0; i < cnt; i ++)
        for (int j = 0; j < sum[i]; j ++)
            ans = mul(ans, primes[i]);
    for (int i = ans.size() - 1; i >= 0; i --) cout << ans[i];
    return 0;
}