python/TangDou/Topic/组合数公式专题.md

C_a^b的多种场景下的求法

一、AcWing 885. 求组合数 I

理论依据：\large C_a^b=C_{a-1}^b+C_{a-1}^{b-1}

限制条件： \large a<=2000,b<=2000

受限原因：需要二维数组，2000*2000=4000000=4e6,占据太大的空间。

证明: 有a个苹果，现在需要选出b个苹果,一共有多少种选法呢？ 答：走到第一个苹果面前，面临两个选择：

• 选择它
• 放弃它

如果选择它，将要面对a-1个苹果中选b-1个苹果的问题

如果放弃它，将要面对a-1个苹果中选b个苹果的问题

根据加法原理，两种选择方法加在一起，就是选法总数，也就是上面的递推式。

#include <bits/stdc++.h>

using namespace std;
const int N = 2010;
const int p = 1e9 + 7;
int c[N][N];

void init() {
    /*
    功能：第0列初始化为1，现实意义：C(n,0)=1
    理解：从n个元素中,每次取出0个元素,也就是每次都不取出元素,有几种取法?
    很明显只有1种取法,就是什么都不取这样一种情况。类似:C(n,0)=C(n,n)=1
    这是递归的base case,如果没有这个1打底，以后再怎么加法原理都是白费
    */
    for (int i = 1; i < N; i++) c[i][0] = 1, c[i][i] = 1; // base case

    for (int i = 1; i < N; i++)
        for (int j = 1; j < i; j++)
            // 递归式最重要！这是根据数学公式得到的递推式
            // 从递推式由右向左看过去，需要由小到大填充i,j,并且，i>j
            c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % p;
}

int n;

int main() {
    init();
    cin >> n;
    while (n--) {
        int a, b;
        cin >> a >> b;
        cout << c[a][b] << endl;
    }
    return 0;
}

二、AcWing 886. 求组合数 II

上个办法把C_a^b的值预处理出来,用的是c[N][N],N最大2010. 本题的a,b都是上限10^5,如果按上题来，就是c[10^5][10^5], 直接报 Memory\ Limit\ Exceeded. 所以不能直接递推求出所有解。

限制条件： a<=1e5,b<=1e5

理论依据：\large C_a^b=\frac{a!}{(a-b)! * b!}

举栗子： \large C_3^2=\frac{3 \times 2}{2 \times 1}=3

前导知识

模逆元

定义与用途

若b,p互质，b \times b^{-1} \equiv 1 \ (mod \ p),称b^{-1}是b的模p逆元。

有 $a \equiv a \ (mod \ p) \ \equiv a \times 1 \ (mod \ p) \equiv a \times (b \times b^{-1}) \ (mod \ p)$

若b \mid a,则\frac{a}{b} \equiv a \times b^{-1} \ (mod \ p)

阶乘模逆元

若

a \times a^{-1} \equiv 1 \ (mod \ p) \\ b \times b^{-1} \equiv 1 \ (mod \ p)

则

ab \times a^{-1}  b^{-1} \equiv 1 \ (mod \ p)

根据逆元定义，则

(ab)^{-1} \equiv a^{-1}b^{-1} \ (mod \ p)

因此

\large \displaystyle (n!)^{-1}= \prod_{i=1}^{n}i^{-1} \ (mod \ p)

本题要求结果 mod\ (1e9+7), 模是质数(一般取模的都是质数)，如果有除法，可以直接使用 费马小定理+快速幂 求逆元转化为乘法求解。

现在要求计算出 a!, {b!}^{-1}, {(a-b)!}^{-1}，可以使用两个数组来递推得到：

\large fact[i] = i!\ \% \ p

\large infact[i] = (i!)^{-1}\ \% \ p = ((i-1)!^{-1}*i^{-1}\ ) \% \ p

解释：上面推导的阶乘模逆元

所以： C_a^b=\frac{a!}{(a-b)! * b!} = (fact[a] * infact[a-b] * infact[b] )\ \% \ p

C++ 代码

#include <bits/stdc++.h>

using namespace std;
#define int long long
#define endl "\n"
const int N = 100010;
const int p = 1e9 + 7;

int fact[N];   // 阶乘结果数组
int infact[N]; // 阶乘逆元结果数组

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

signed main() {
    // base case
    fact[0] = infact[0] = 1;

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

    int n;
    cin >> n;
    while (n--) {
        int a, b;
        cin >> a >> b;
        printf("%lld\n", fact[a] * infact[b] % p * infact[a - b] % p);
    }
}

AcWing 887. 求组合数 III

数据范围：b<=a<=1e18,p<=1e5,p是质数

本题中a,b都太大了，如果像上一题一样预处理出fact[],infact[],直接就炸了，肯定行不通！

还好，它的p比较小，而且是质数，我们有一种武器叫Lucas公式：

1.Lucas公式

\large C_a^b \equiv C_{a\%p}^{b\%p} * C_{a \div p}^{b \div p} (mod \ p)

注意： 这个p是小于1e5的，看到要求模1e9+7可不能用Lucas定理啊！不是干那个的啊！适合a,b很大，p很小的场景啊！！

2.此种情况下如何求解C_a^b

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

根据这个变形:

• 对于a来讲，需要变量j从a一直遍历到a-b+1
• 对于b来讲，需要变量i从1一起遍历到b

从a倒着走到a-b+1是b次(a>=b)，比如10遍历到8，就是10,9,8三次，即b=3,也就是10-3+1=8。 而i也是遍历b次，太巧了！它们两个可以在一个b次的循环中一起变动 (代码能省则省啊~)！

int C(int a, int b) {
    if (a < b) return 0;
    int down = 1, up = 1;
    for (int i = a, j = 1; j <= b; i--, j++) {
        up = up * i % p;
        down = down * j % p;
    }
    return up * qmi(down, p - 2, p) % p;
}

Code

#include <bits/stdc++.h>

using namespace std;
#define int long long
#define endl "\n"

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

//  功能：组合数模板
int C(int a, int b, int p) {
    if (a < b) return 0;
    int down = 1, up = 1;

    for (int i = a, j = 1; j <= b; i--, j++) {
        up = up * i % p;
        down = down * j % p;
    }
    return up * qmi(down, p - 2, p) % p;
}

//  Lucas公式
int lucas(int a, int b, int p) {
    if (a < p && b < p) return C(a, b, p);
    return C(a % p, b % p, p) * lucas(a / p, b / p, p) % p; // 递归套用公式
}

int n, p;

signed main() {
    cin >> n;
    while (n--) { // n组询问
        int a, b;
        cin >> a >> b >> p;
        cout << lucas(a, b, p) << endl; // 利用lucas公式计算组合数
    }
    return 0;
}

四、AcWing 888. 求组合数 IV

不让取模怎么办?

1、与前三道题的区别

不再是数据上限的问题了，而是一直不mod,有多大都保留，这无疑可能会爆long\ long,需要使用高精度。

2、公式一【从定义出发的组合数公式】：

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

3、是不是我们利用高精度+上面的组合数公式直接算就行了？

不是的，因为yxc(大雪菜老师)说，这样的效率太低，不可以，需要再想一个好办法。

4、公式二 【算术基本定理】

算术基本定理:\large C_a^b=p_1^{\alpha_1} \times p_2^{\alpha_2} \times p_3^{\alpha_3} ... \times p_k^{\alpha_k} 其中p_1,p_2...是质数，\alpha 1,\alpha 2,...是指质数p_1,p_2,...的个数。 如果能分解质因数成功的话，那么就可以通过 高精度乘法 解决掉这个问题。

(1) 那么p_1和p_2这些东东怎么求？ 思路：因为a,b都是在[1,5000]之间的数，所以可以通过筛质数的方法，提前获取到一个质数数组，然后逐个去看，是不是含有这个质数。就能知道有哪些p_1,p_2,...了。

(2) \alpha_1,\alpha_2,...,\alpha_k又该怎么求呢？

5、公式三【求质数p在a!中出现的次数】

\large cnt=\lfloor \frac{a}{p} \rfloor + \lfloor \frac{a}{p^2} \rfloor + \lfloor \frac{a}{p^3} \rfloor + ...+ \lfloor \frac{a}{p^k} \rfloor

参考资料

计算n的阶乘中质数p的个数：

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

举个栗子：

5!=1 \times 2 \times 3 \times 4 \times 5,求2的个数， 则1 \sim 5之间含2的个数就是 \lfloor \frac{5}{2} \rfloor,就是2个。 则1 \sim 5之间含2^2的个数就是 \lfloor \frac{5}{2^2} \rfloor,就是1个。 2^3就大于5，就不再继续。那么一共的个数就是2+1=3个。表示在 5!这个数字中，存在3个2，就是2中有1个2，4中有两个2。

6、算法流程

(1)、筛素数，把1-5000之间的素数筛出来。 (2)、计算C_a^b中每个已求得素数的个数。 (3)、利用高精度，计算C_a^b=p_1^{\alpha1} \times p_2^{\alpha2} \times p_3^{\alpha3} ... \times p_k^{\alpha k}的值。

Code

#include <bits/stdc++.h>

using namespace std;
const int N = 5010;

int primes[N], cnt;
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;
        }
    }
}

// 高精乘低精
void mul(int a[], int &al, int b) {
    int t = 0;
    for (int i = 1; i <= al; i++) {
        t += a[i] * b;
        a[i] = t % 10;
        t /= 10;
    }
    while (t) {
        a[++al] = t % 10;
        t /= 10;
    }
}

/**
 * 功能：n的阶乘中包含的质因子p的个数
 * @param n n的阶乘
 * @param p 质因子p
 * @return 有多少个
 */
int get(int n, int p) {
    int cnt = 0;
    while (n) { // p^1的个数，p^2的个数,p^3的个数...
        cnt += n / p;
        n /= p;
    }
    return cnt;
}

// C(a,b)的结果，高精度保存到c数组，同时，返回c数组的长度len
void C(int a, int b, int c[], int &cl) {
    // 乘法的基数是1
    c[1] = 1, cl = 1; // 由于高精度数组中只有一位，是1，所以长度也是1

    for (int i = 0; i < cnt; i++) { // 枚举区间内所有质数
        int p = primes[i];
        /*
        C(a,b)=a!/(b! * (a-b)!)
        a!中有多少个质数因子p
        减去(a-b)!的多少个质数因子p,
        再减去b!的质数因子p的个数，就是总个数
        s记录了p这个质数因子出现的次数
        */
        int s = get(a, p) - get(b, p) - get(a - b, p);
        while (s--) mul(c, cl, p); // 不断的乘p,结果保存到数组c中。len将带回c的有效长度
    }
}

int a, b;
int c[N], cl;

int main() {
    cin >> a >> b;
    // 筛质数
    get_primes(N - 1);

    // 计算组合数
    C(a, b, c, cl);

    // 输出高精度结果
    for (int i = cl; i >= 1; i--) printf("%d", c[i]);
    return 0;
}

7、这段代码如何理解？

 s = get(a, p) - get(a - b, p) - get(b, p); 

看公式一： C_a^b=\frac{a \times (a-1) \times ... \times(a-b+1)}{b\times (b-1) \times ... \times 1} =\frac{a \times (a-1) \times ... \times (a-b+1) \times (a-b)!}{ b \times (b-1) \times ... \times 1 \times (a-b)!} =\frac{a!}{b! \times (a-b)!} a!中质数p的个数，减去b!中质数p的个数，再减去(a-b)!中质数p的个数，就是公因子消掉的意思。