组合数学 —— 组合数取模 —— 卢卡斯定理与扩展卢卡斯定理

【卢卡斯定理】

1.要求：p 是质数，m、n 很大但 p 很小 或者 n、m 不大但大于 p

2.定理内容

• 
•  

其中，

3.推论

当将 n 写成 p 进制：，将 m 写成 p 进制： 时，有：

• 

4.实现

代码实现可简单理解为：

LL fac[N];
void getFac(){//构造阶乘
    fac[0]=1;
    for(int i=1;i<1000000;i++){
        fac[i]=fac[i-1]*i%MOD;
    }
}
LL quickPowMod(LL a,LL b,LL mod){//快速幂
    LL res=1;
    while(b){
        if(b&1)
            res=res*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return res;
}
LL getC(LL n,LL m,LL mod){//获取C(n,m)%mod
    if(m>n)
        return 0;
    return fac[n]*(quickPowMod(fac[m]*fac[n-m]%mod,mod-2,mod))%mod;
}
LL Lucas(LL n,LL m,LL mod){//卢卡斯定理
    if(m==0)
        return 1;
    return getc(n%mod,m%mod,mod)*Lucas(n/mod,m/mod,mod)%mod;
}
int main(){
    getFac();
    LL n,m;
    scanf("%lld%lld",&n,&m);
    printf("%lld\n",Lucas(n,k,MOD));
    return 0;
}

【扩展卢卡斯】

卢卡斯定理适用于 p 是素数的情况，但当 p 不是素数时，可以将其分解质因数，将组合数按照卢卡斯定理的方法求 p 的质因数的模，然后用中国剩余定理合并即可。

要求：p 不是质数 且 m、n 很大但 p 很小 或者 n、m 不大但大于 p

例如： 

当需要计算   时，可以求出：

然后对于方程组：，可以求出满足条件的最小的 x，记为：

那么有：

但是， 并不是一个素数，而是某个素数的某次方，那么就需要计算  
对于 ，已知 ，因此若能计算出 ，就能计算出  和 

设 ，那么答案就是 ，其中  代表计算 a 对 b 的乘法逆元

于是，问题就转换为如何计算 

例如：p=3，t=2，n=19，有：

n! = 1×2×3×4×5×6×7×8×…×19

    = (1×2×4×5×7×8×…×16×17×19) × (3×6×9×12×15×18)

    = (1×2×4×5×7×8×…×16×17×19) × 36 × (1×2×3×4×5×6) 

后半部分是 ，递归即可。前半部分是以  为周期的

下面是孤立的 19，可以知道孤立出来的长度不超过 ，直接计算即可。

对于最后剩下的 36，只要计算出 n!、m!、(n−m)! 里含有多少个 p ，设他们分别有 x、y、z 个 p，那么 x−y−z 就是   中 p 的个数，直接计算即可

LL powerMod(LL a,LL b,LL p) {//快速幂取模
    LL ans=1;
    while(b){
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans;
}
LL fac(LL n,LL p,LL pk) {//计算阶乘
    if(!n)
        return 1;
    LL ans=1;
    for(int i=1; i<pk; i++)
        if(i%p)
            ans=ans*i%pk;
    ans=powerMod(ans,n/pk,pk);
    for(int i=1;i<=n%pk;i++)
        if(i%p)
            ans=ans*i%pk;
    return ans*fac(n/p,p,pk)%pk;
}
LL extendGCD(LL a,LL b,LL &x,LL &y) {//扩展欧几里得
    if(!b){
        x=1;
        y=0;
        return a;
    }
    LL xx, yy;
    LL gcd=extendGCD(b,a%b,xx,yy);
    x=yy;
    y=xx-a/b*yy;
    return gcd;
}
LL inv(LL a,LL p){//计算逆元
    LL x,y;
    extendGCD(a,p,x,y);
    return (x%p+p)%p;
}
LL C(LL n,LL m,LL p,LL pk) {//组合数模质数幂
    if(n<m)
        return 0;
    LL f1=fac(n,p,pk);
    LL f2=fac(m,p,pk);
    LL f3=fac(n-m,p,pk);
    LL cnt=0;
    for(LL i=n; i; i/=p)
        cnt+=i/p;
    for(LL i=m; i; i/=p)
        cnt-=i/p;
    for(LL i=n-m; i; i/=p)
        cnt-=i/p;
    return f1 * inv(f2,pk)%pk * inv(f3, pk)%pk * powerMod(p,cnt,pk)%pk;
}
LL a[N],c[N];
int tot;
LL CRT() {//中国剩余定理
    LL M=1,ans=0;
    for (int i=0; i<tot; i++)
        M*=c[i];
    for (int i=0; i<tot; i++)
        ans=(ans + a[i] * (M/c[i])%M * inv(M/c[i],c[i])%M )%M;
    return ans;
}
LL extendLucas(LL n,LL m,LL p) {//扩展卢卡斯
    for(int i=2; p>1&&i<=sqrt(p); i++) {
        LL temp=1;
        while(p%i==0) {
            p/=i;
            temp*=i;
        }
        if(temp>1) {
            a[tot]=C(n,m,i,temp);
            c[tot++]=temp;
        }
    }
    if(p>1) {
        a[tot]=C(n,m,p,p);
        c[tot++]=p;
    }
    return CRT();
}
int main() {
    LL n,m,p;
    cin>>n>>m>>p;
    cout<<extendLucas(n,m,p);
    return 0;
}