Lucas 定理

简述

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

$Lucas$

$Lucas$ 定理内容如下：
对于质数 $P$，有

$$(\genfrac{}{}{0ex}{}{n}{m})\equiv (\genfrac{}{}{0ex}{}{\lfloor \frac{n}{P}\rfloor }{\lfloor \frac{m}{P}\rfloor })\cdot (\genfrac{}{}{0ex}{}{n\mathrm{mod}P}{m\mathrm{mod}P})\mathrm{mod}P$$

上式中 $n\mathrm{mod}P$ 和 $m\mathrm{mod}P$ 一定小于$P$，可以直接求解，$(\genfrac{}{}{0ex}{}{\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}$$

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

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

$$\{\begin{array}{rl}{a}_1\equiv {\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}& (\mathrm{mod}{p_1}^{{\alpha }_1})\\ a_2\equiv {\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}& (\mathrm{mod}{p_2}^{{\alpha }_2})\\ & \cdots \\ a_r\equiv {\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}& (\mathrm{mod}{p_r}^{{\alpha }_r})\end{array}$$

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

根据同余的定义，${\displaystyle a_i=(\genfrac{}{}{0ex}{}{n}{m})mod{p_i}^{{\alpha }_i}}$，问题转化成，求 ${\displaystyle (\genfrac{}{}{0ex}{}{n}{m})mod{p}^{\alpha }}$（$p$ 为质数）的值。

根据组合数定义 ${\displaystyle (\genfrac{}{}{0ex}{}{n}{m})=\frac{n!}{m!(n-m)!}}$，${\displaystyle (\genfrac{}{}{0ex}{}{n}{m})mod{p}^{\alpha }=\frac{n!}{m!(n-m)!}mod{p}^{\alpha }}$。

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

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

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

所以将原式转化为

$$\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}{p}^{x-y-z}mod{p}^{\alpha }$$

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

问题转化成求

$$\frac{n!}{q^x}mod{q}^k$$

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

得到

$$n!=q^{\lfloor \frac{n}{q}\rfloor }\cdot \left(\lfloor \frac{n}{q}\rfloor \right)!\cdot {\left(\prod _{i,(i,q)=1}^{q^k}i\right)}^{\lfloor \frac{n}{q^k}\rfloor }\cdot \left(\prod _{i,(i,q)=1}^{nmod{q}^k}i\right)$$

于是

$$\frac{n!}{q^{\lfloor \frac{n}{q}\rfloor }}=\left(\lfloor \frac{n}{q}\rfloor \right)!\cdot {\left(\prod _{i,(i,q)=1}^{q^k}i\right)}^{\lfloor \frac{n}{q^k}\rfloor }\cdot \left(\prod _{i,(i,q)=1}^{nmod{q}^k}i\right)$$

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

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

$$\frac{n!}{q^x}=\frac{\left(\lfloor \frac{n}{q}\rfloor \right)!}{q^{x^{\prime }}}\cdot {\left(\prod _{i,(i,q)=1}^{q^k}i\right)}^{\lfloor \frac{n}{q^k}\rfloor }\cdot \left(\prod _{i,(i,q)=1}^{nmod{q}^k}i\right)$$

其中 $x$ 和 $x^{\prime }$ 都表示把分子中所有的素数 $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