$Lucas$ 定理

简述

当问题规模很大，而模数是一个不大的质数的时候，就不能简单地通过递推求解来得到答案，需要用到 Lucas 定理。

\(Lucas\)

\(Lucas\) 定理内容如下：
对于质数 \(P\)，有

\[\binom{n}{m} \equiv \binom{\lfloor\frac{n}{P}\rfloor}{\lfloor\frac{m}{P}\rfloor}\cdot\binom{n\mod P}{m\mod P}\mod P \]

上式中 \(n\mod P\) 和 \(m\mod P\) 一定小于\(P\)，可以直接求解，\(\binom{\lfloor\frac{n}{P}\rfloor}{\lfloor\frac{m}{P}\rfloor}\)继续用 \(Lucas\) 定理求解。边界条件：当 \(m=0\) 的时候，返回 \(1\)
复杂度为 \(O(f(p) + g(n)\log n)\)，其中 \(f(n)\) 为预处理组合数的复杂度，\(g(n)\) 为单次求组合数的复杂度。

il int Lucas(int n,int m,int p){
    if(!m) return 1;
    return 1ll*Lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}

code

#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
using namespace std;
il int rd(){
    ri int x=0,f=1;ri char c=getchar();
    while(c<'0'||c>'9') f=(c=='-')?-1:1,c=getchar();
    while(c>='0'&&c<='9') x=x*10+(c^48),c=getchar();
    return x*f;
}
cs int N=1e5+5;
int t,n,m,P,A[N];
il int qpow(int a,int b,int p){
    ri int ans=1;
    while(b) {
        if(b&1) ans=(1ll*ans*a)%p;
        a=(1ll*a*a)%p,b>>=1;
    }
    return ans;
}
il int inv(int a,int p){
    return qpow(a,p-2,p);
}
il int C(int n,int m,int p){ // 直接算组合数
    if(m>n) return 0;
    return 1ll*A[n]*inv(A[m],p)%p*inv(A[n-m],p)%p;
}
il int Lucas(int n,int m,int p){
    if(!m) return 1;
    return 1ll*Lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
int main(){
    t=rd();
    while(t--){
        n=rd(),m=rd(),P=rd(),A[1]=A[0]=1; // 预处理阶乘
        for(ri int i=2;i<P;++i) A[i]=1ll*A[i-1]*i%P;
        printf("%d\n",Lucas(n+m,n,P));
    }
    return 0;
} 

\(exLucas\)

对于 \(P\) 不是素数的情况，就需要用到 \(exLucas\)

将 \(P\) 质因数分解：

\[P={p_1}^{\alpha_1}\cdot{p_2}^{\alpha_2}\cdots{p_r}^{\alpha_r}=\prod_{i=1}^{r}{p_i}^{\alpha_i} \]

只需求出 \(\binom{n}{m}\bmod {p_i}^{\alpha_i}\) 的值（其中 \(p_i\) 为素数且 \(\alpha_i\) 为正整数），利用 中国剩余定理 合并答案即可。

具体的，对于任意 \(i,j\)，有 \({p_i}^{\alpha_i}\) 与 \({p_j}^{\alpha_j}\) 互质，可以构造如下 \(r\) 个同余方程

\[\left\{ \begin{aligned} a_1\equiv \displaystyle\binom{n}{m}&\pmod {{p_1}^{\alpha_1}}\\ a_2\equiv \displaystyle\binom{n}{m}&\pmod {{p_2}^{\alpha_2}}\\ &\cdots\\ a_r\equiv \displaystyle\binom{n}{m}&\pmod {{p_r}^{\alpha_r}}\\ \end{aligned} \right. \]

于是在求出 \(a_i\) 后，就可以用中国剩余定理求解出 \(\displaystyle\binom{n}{m}\)。

根据同余的定义，\(\displaystyle a_i=\binom{n}{m}\bmod {p_i}^{\alpha_i}\)，问题转化成，求 \(\displaystyle \binom{n}{m} \bmod p^\alpha\)（\(p\) 为质数）的值。

根据组合数定义 \(\displaystyle \binom{n}{m} = \frac{n!}{m! (n-m)!}\)，\(\displaystyle \binom{n}{m} \bmod p^\alpha = \frac{n!}{m! (n-m)!} \bmod p^\alpha\)。

由于式子是在模 \(p^\alpha\) 意义下，所以分母要算乘法逆元。

同余方程 \(ax \equiv 1 \pmod p\) 有解的充要条件为 \(\gcd(a,p)=1\)（裴蜀定理），

然而无法保证有解，发现无法直接求 \(\operatorname{inv}_{m!}\) 和 \(\operatorname{inv}_{(n-m)!}\)，

所以将原式转化为

\[\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}p^{x-y-z} \bmod p^\alpha \]

\(x\) 表示 \(n!\) 中包含多少个 \(p\) 因子，\(y\), \(z\) 同理。

问题转化成求

\[\frac{n!}{q^x}\bmod q^k \]

这时可以利用 Wilson定理 的推论。

得到

\[n! = q^{\left\lfloor\frac{n}{q}\right\rfloor} \cdot \left(\left\lfloor\frac{n}{q}\right\rfloor\right)! \cdot {\left(\prod_{i,(i,q)=1}^{q^k}i\right)}^{\left\lfloor\frac{n}{q^k}\right\rfloor} \cdot \left(\prod_{i,(i,q)=1}^{n\bmod q^k}i\right) \]

于是

\[\frac{n!}{q^{\left\lfloor\frac{n}{q}\right\rfloor}} = \left(\left\lfloor\frac{n}{q}\right\rfloor\right)! \cdot {\left(\prod_{i,(i,q)=1}^{q^k}i\right)}^{\left\lfloor\frac{n}{q^k}\right\rfloor} \cdot \left(\prod_{i,(i,q)=1}^{n\bmod q^k}i\right) \]

\(\displaystyle \left(\left\lfloor\frac{n}{q}\right\rfloor\right)!\) 同样是一个数的阶乘，所以也可以分为上述三个部分，于是可以递归求解。

等式的右边两项不含素数 \(q\)。如果直接把 \(n\) 的阶乘中所有 \(q\) 的幂都拿出来，等式右边的阶乘也照做，该式可以直接写成

\[\frac{n!}{q^{x}} = \frac{\left(\left\lfloor\frac{n}{q}\right\rfloor\right)!}{q^{x'}} \cdot {\left(\prod_{i,(i,q)=1}^{q^k}i\right)}^{\left\lfloor\frac{n}{q^k}\right\rfloor} \cdot \left(\prod_{i,(i,q)=1}^{n\bmod q^k}i\right) \]

其中 \(x\) 和 \(x'\) 都表示把分子中所有的素数 \(q\) 都拿出来。这样，每一项就完全不含 \(q\) 了。

时间复杂度 \(O(p\log p)\)

code

#include<bits/stdc++.h>
#define il inline int
#define vl inline void
#define cs const
#define ri register
#define int long long
using namespace std;
vl exgcd(int a,int b,int &x,int &y){
    (b==0)?(x=1,y=0):(exgcd(b,a%b,y,x),y-=(a/b)*x);
}
il inv(int a,int p){
    int x,y;exgcd(a,p,x,y);
    return (x%p+p)%p;
}
il qpow(int a,int b,int p){
    int as=1;
    while(b){
        if(b&1) as=(as*a)%p;
        a=(a*a)%p,b>>=1;
    }
    return as;
}
il getf(int n,int p,int pk){
    if(!n) return 1;
    int as=1;
    for(ri int i=2;i<pk;++i){
        if(i%p) as=(as*i)%pk;
    }
    as=qpow(as,n/pk,pk);
    for(ri int i=1+pk*(n/pk);i<=n;++i){
        if(i%p) as=(i%pk*as)%pk;
    }
    return as*getf(n/p,p,pk)%pk;
}
il getg(int n,int p){
    return (n<p)?0:(getg(n/p,p)+n/p);
}
il getc(int n,int m,int p,int pk){
    int fz=getf(n,p,pk),fm1=getf(m,p,pk),fm2=getf(n-m,p,pk);
    int qp=qpow(p,getg(n,p)-getg(m,p)-getg(n-m,p),pk);
    return qp*fz%pk*inv(fm1,pk)%pk*inv(fm2,pk)%pk;
}
il exlucas(int n,int m,int p){
    int res=p,as=0,pk;
    for(ri int i=2;i*i<p;++i){
        if(res%i==0){
            pk=1;
            while(res%i==0){
                res/=i,pk*=i;
            }
            as=(as+getc(n,m,i,pk)*inv(p/pk,pk)%p*(p/pk)%p)%p;
        }
    }
    if(res>1){
        as=(as+getc(n,m,res,res)*inv(p/res,res)%p*(p/res)%p)%p;
    }
    return as;
}
signed main(){
    int n,m,p;
    cin>>n>>m>>p;
    cout<<exlucas(n,m,p);
    return 0;
} 

edit