schtonn

schtonn

0%

求组合数

Posted on Edited on Valine:

前置知识

中国剩余定理，

引入

从 \(n\) 个不同元素中取出 \(m\) 个元素的组合数，即 \(C_n^m\)，是从 \(n\) 个不同元素中，取出 \(m\) 个元素所形成的组合个数。

我们需要使用各种方式求 \(C_n^m\)，每个方式难度，速度和空间都有不同。

朴素算法

首先，求组合数有一个基本公式：

\[\large C_n^m=\dfrac{P_n^m}{P_m}=\dfrac{n!}{m!(n-m)!}\]

\[\large C_n^0=1\]

我们的朴素算法，就是每遇到一个组合数，都直接套用这个公式求阶乘。

1
2
3
4
5
6
7
8
9
10
11
int fact(int x){
    int ans=1;
    for(int i=2;i<=x;i++){
        ans*=i;
    }
    return ans;
}
int C(int n,int m){
    if(m==0)return 1;
    return fact(n)/(fact(m)*fact(n-m));
}

这个算法适合 \(n,m\leq15\)，时间复杂度极高，但是实现难度为零。

杨辉三角

如果把组合数转换成实际问题，可以考虑每一步选定一个元素，那么一个组合数一定可以由更小的组合数转移而来。

具体思路为，\(C_n^m\) 是 \(n\) 中选 \(m\) 个，考虑从 \(n-1\) 个中选的所有情况，一部分是选新加入的这第 \(n\) 个，一部分是不选。

那么可以列出这样一个式子：

\[\large C_n^m=C_{n-1}^{m-1}+C_{n-1}^m\]

（当然，通过数学的方式推也能推出来）

这时我们发现，如果把这个式子列成表，竖着拎起来，那就正好是一个杨辉三角（每一位等于上面左右两位的和）。

那么根据这个式子，就可以写出一种动态规划算法。

1
2
3
4
5
6
7
8
9
10
11
12
int C[N][M];
void calc(int n,int m){
    for(int i=0;i<=n;i++){
        f[i][i]=1;
        f[i][0]=1;
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<i;j++){
            f[i][j]=f[i-1][j]+f[i-1][j-1];
        }
    }
}

这个算法适合 \(n,m\leq10000\)，时间复杂度适中，需要一个数组，而且实现难度也很低。

Lucas 定理

求 \(C_m^n\bmod p\)，\(p\) 是质数。

对质数 \(p\)，有 Lucas 定理： \[\large C_n^m\bmod p=C_{\lfloor n/p\rfloor}^{\lfloor m/p\rfloor}\cdot C_{n\bmod p}^{m\bmod p}\bmod p\]

证明

（摘自 OI-wiki）

考虑 \(C_p^n\bmod p\) 的取值，注意到 \(C_p^n=\dfrac{p!}{n!(p-n)!}\)，分子的质因子分解中 \(p\) 次项恰为 \(1\)，因此只有当 \(n=0\) 或 \(n=p\) 的时候 \(n!(p-n)!\) 的质因子分解中含有 \(p\)。进而我们可以得出： \[\begin{aligned}(a+b)^p&=\sum_{n=0}^p C_p^n a^nb^{p-n}\\ &\equiv\sum_{n=0}^p[n=0\vee n=p]a^nb^{p-n}\\ &\equiv a^p+b^p\pmod p\end{aligned}\]

注意过程中没有用到费马小定理，因此这一推导不仅适用于整数，亦适用于多项式。因此我们可以考虑二项式 \(f(x)=(ax^n+bx^m)^p\bmod p\) 的结果： \[\begin{aligned}(ax^n+bx^m)^p&\equiv a^px^{pn}+b^px^{pm}\\ &\equiv ax^{pn}+bx^{pm}\\ &\equiv f(x^p)\end{aligned}\]

考虑二项式 \((1+x)^n\mod p\)，那么 \(C_n^m\) 就是求其在 \(x^m\) 次项的取值。使用上述引理，我们可以得到： \[\begin{aligned}(1+x)^n&\equiv(1+x)^{p\lfloor n/p\rfloor}(1+x)^{n\bmod p}\\&\equiv(1+x^p)^{\lfloor n/p\rfloor}(1+x)^{n\bmod p}\end{aligned}\]

注意前者只有在 \(p\) 的倍数位置才有取值，而后者最高次项为 \(n\bmod p\leq p-1\)，因此这两部分的卷积在任何一个位置只有最多一种方式贡献取值，即在前者部分取 \(p\) 的倍数次项，后者部分取剩余项，即 \(C_n^m\bmod p=C_{\lfloor n/p\rfloor}^{\lfloor m/p\rfloor}\cdot C_{n\bmod p}^{m\bmod p}\bmod p\)。

求解

观察 Lucas 定理，可知 \(n\bmod p\) 和 \(m\bmod p\) 一定是小于 \(p\) 的数，可以直接求解 \(C_{\lfloor n/p\rfloor}^{\lfloor m/p\rfloor}\)，递归使用 Lucas 定理。

\(C_{n\bmod p}^{m\bmod p}\) 则由于 \(p\) 是质数，利用费马小定理转化成 \(n!\) 乘 \(m!(n-m)!\) 的逆元，也就是 \(m!(n-m)!^{p-2}\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
long long pow(long long a,long long b,long long m){
    long long ans=1;
    a%=m;
    while(b){
        if(b&1)ans=(ans%m)*(a%m)%m;
        b/=2;
        a=(a%m)*(a%m)%m;
    }
    ans%=m;
    return ans;
}
long long inv(long long x,long long p){
    return pow(x,p-2,p);
}
long long C(long long n,long long m,long long p){
    if(m>n)return 0;
    long long up=1,down=1;
    for(int i=n-m+1;i<=n;i++)up=up*i%p;
    for(int i=1;i<=m;i++)down=down*i%p;
    return up*inv(down,p)%p;
}
long long Lucas(long long n,long long m,long long p){
  if(m==0)return 1;
  return(Lucas(n/p,m/p,p)*C(n%p,m%p,p))%p;
}

这个算法适合 \(n,m\leq10^8,p\leq10^4\)，速度较快，理解后难度也不高。

拓展 Lucas 定理

适用于求 \(C_m^n\bmod p\)，\(p\) 不一定是质数。 ### 分解 将 \(p\) 分解质因数： \[p=q_1^{a_1}\cdot q_2^{a_2}\cdots q_r^{a_r}=\prod_{i-1}^rq_i^{a_i}\]

对于任意 \(i,j\)，有 \(p_i^{a_i}\) 与 \(p_j^{a_j}\) 互质，所以可以构造 \(r\) 个同余方程：

\[ \left\{ \begin{aligned} a_1&\equiv C_n^m\pmod {q_1^{a_1}}\\ a_2&\equiv C_n^m\pmod {q_2^{a_2}}\\ &\cdots\\ a_r&\equiv C_n^m\pmod {q_r^{a_r}} \end{aligned} \right. \]

我们发现，在求出 \(a_i\) 后，就可以用中国剩余定理求解出 \(C_n^m\)。

转化

根据同余的定义，\(a_i=C_n^m\bmod q_i^{a_i}\)，问题转化成，求 \(C_n^m\bmod\)