求组合数（四种方法）

求组合数（四种方法）

文章首发于我的个人博客：欢迎大佬们来逛逛

递推（杨辉三角）

对于求一个数字的组合数：$C_n^m$ 可以分解为这两种形式：

• $C_{n-1}^{m-1}$ ：选择当前数
• $C_{n-1}^m$ ：不选择当前数

因此就得到了递推公式：

$$C_n^m=C_{n-1}^{m-1}+C_{n-1}^m$$

for (int i=0;i<=n;i++){
        for (int j=0;j<=i;j++){
            if (j==0){
                dp[i][j]=1; //n个数中取0个，默认为1
            }
            else{
                dp[i][j]=(dp[i-1][j]+dp[i-1][j-1])%p;
            }
        }
    }

由于需要用到二重循环，因此此方法适用的范围是： $n,m<=10^4$ 左右

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

---

快速幂+乘法逆元

考虑组合数的计算公式：

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

我们可以直接计算这个式子中的每一项。

费马小定理：如果 $p$ 是一个质数，并且 $a$ 和 $p$ 互质，则满足： $a^{p-1}\equiv 1(modp)$

乘法逆元：如果 $a$ 与 $p$ 都是质数，并且满足同余方程式：$a\cdot x\equiv 1(modp)$ ，则称 $x$ 为 $a$ 模 $p$ 的乘法逆元，即：$\frac{1}{amodp}=x$

---

那么我们可以将费马小定理进行转化，得到： $a\cdot a^{p-2}\equiv 1(modp)$

则如果 $a$ 和 $p$ 都是质数，$x=a^{p-2}$ 就是 $a$ 模 $p$ 的乘法逆元。

因此我们要求 $a$ 模 $p$ 的乘法逆元，直接计算 $a^{p-2}$ 即可

---

我们规定：

• $f$ 数组为阶乘数组，保存每个数字阶乘意义下的阶乘。
• $g$ 数组为逆元数组，保存每个数组阶乘意义下的逆元。

因此就可以转换为 $C_n^m=f[n]\cdot g[n-m]\cdot g[m]$

如何求这个数组的值呢？

• 对于阶乘来说比较容易。
• 对于逆元，需要进行转化： $\frac{1}{i!}=\frac{1}{(i-1)!}\cdot \frac{1}{i}$ ，因此就是$g[i]=g[i-1]\cdot inv(i)$

f[0]=g[0]=1; //初始值
for (int i=1;i<=n;i++){
    f[i]=f[i-1]*i%p; //计算i的阶乘
    g[i]=g[i-1]*qpow(i,p-2)%p; //计算i的乘法逆元 qpow为快速幂
}
...
int get(int n,int m){ //得到C(n,m)的组合数答案
    return f[n]*g[m]*g[n-m]%p;
}

容易得知此方法计算组合数的使用范围大概是： $n,m<=10^7$ ，其主要限制为 $for$ 循环和数组空间大小

时间复杂度：$O(nlo{g}_{2}n)$

---

卢卡斯定理

此方法用于非常大的 $n$ 和 $m$ ，但是 $p$ 却比较小的计算。

定理：

$$C_n^m=C_{\frac{n}{p}}^{\frac{m}{p}}\cdot C_{n\%p}^{m\%p}\%p$$

即计算组合数就可以利用拆分的方法。（就不证明了）

然后注意到 $p$ 是一个质数，因此 $n$ 和 $m$ 对其取余后，最后的结果一定比 $p$ 小。

因此就可以对第二个式子进行 快速幂+逆元 的方法，我们知道这个方法的瓶颈是过大的 $n$ 导致空间不够或者循环过大，通过这种方法缩小后便可以直接使用此方法。

我们便可以利用递归来依次缩小 $n$ 和 $m$ 的值，进而计算部分组合数，最后相乘即可，递归结束的条件是 $m==0$ ，此时返回 1 即可。

注意在计算某个 $C_i^j$ 的时候，我们的两个数组的预处理是到 $p$ 为止（因为不可能比 $p$ 还大了）的。

f[0]=g[0]=1; //初始值
for (int i=1;i<=p;i++){
    f[i]=f[i-1]*i%p; //计算i的阶乘
    g[i]=g[i-1]*qpow(i,p-2)%p; //计算i的乘法逆元 qpow为快速幂
}
...
int get(int n,int m){ //得到C(n,m)的组合数答案
    return f[n]*g[m]*g[n-m]%p;
}
...
int lucas(int n,int m){
    if (m==0){
        return 1;
    }
    return lucas(n/p,m/p)*get(n%p,m%p)%p;
}

此方法使用范围（大概）： $n,m<=10^{18}，p<=10^6$

时间复杂度：$O(n\cdot lo{g}_{2}n)$

---

高精度组合数

如果我们需要求 $C_n^m$ 而不需要求其模数，我们可以发现这个数字将会变得非常大，因为组合数的增长是非常快的！此时我们必须要采用高精度的方法。

此方法我们仍然需要用到公式：

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

步骤：

1. 首先预处理所有 $[1,n]$ 范围内的质数，使用线氏筛。
2. 然后统计某个数阶乘意义下能够包含质数的个数
   1. 先处理 $n!$ 能够包含某质数的个数
   2. 再处理 $(n-m)!$ 能够包含某质数的个数
   3. 再处理 $m!$ 能够包含某质数的个数
   4. 最后一式减去二三两式就处理出了 $C_i^j$ 所包含的某质数的个数
3. 最后对所有的质数都执行相同的操作，然后用 $C_i^j\cdot p$ ，其中 p是质数。

其中 $C_i^j\cdot p$ 可以表示为：

$$\prod _a^{b}C[i]\cdot p,p\in [1,包含的p的个数]$$

我们用高精度乘法来维护。

int get(int n,int p){
    //n!有多少p质数
    int cnt=0;
    while (n){
        cnt+=n/p;
        n/=p;
    }
    return cnt;
}
int how(int n,int m,int p){
    return get(n,p)-get(n-m,p)-get(m,p);
}
void mul(int s,int& len){
    //高精度乘法
    int t=0;
    for(int i=0;i<len;i++){
        t+=C[i]*s;
        C[i]=t%10;
        t/=10;
    }
    while (t){
        C[len++]=t%10;
        t/=10;
    }
}
signed main(){
    //求C(n,m)
    std::cin>>n>>m; 
    init();
    int len=1;
    C[0]=1;
    for (int i=1;i<=cnt;i++){
        int x=primes[i];
        int s=how(n,m,x); //C[i]中有多少质数x
        while (s--){
            mul(x,len);
        }
    }
    for (int i=0;i<len;i++){
        std::cout<<C[len-i-1];
    }
    return 0;
}

---