卢卡斯定理

Lucas定理

用于求解大组合数取模且模数为质数的情况。

定理，\(\binom{n}{m}=\binom{\left \lfloor \frac{n}{p}\right \rfloor}{\left \lfloor \frac{m}{p}\right \rfloor}*\binom{n\mod p}{m\mod p}\).

边界条件，当 \(m=0\) 时返回 \(1\).

一般模数的范围不超过 \(10^5\).

P3807 【模板】卢卡斯定理/Lucas 定理

每次给定 \(n,m,p,\)，求 \(\binom{n+m}{m}(\mod p)\).

预处理阶乘到 \(p-1\)，之后进行卢卡斯分解求组合数。

int Lucas(int x,int y){
    if(!y)return 1;
    return Lucas(x/mod,y/mod)*C(x%mod,y%mod)%mod;
}
signed main(){
    int t;
    cin>>t;
    while(t--){
        cin>>n>>m>>mod;
        fac[0]=1;
        for(int i=1;i<=mod-1;i++)fac[i]=fac[i-1]*i%mod;
        facinv[mod-1]=po(fac[mod-1],mod-2,mod);
        for(int i=mod-1;i;i--)facinv[i-1]=facinv[i]*i%mod;
        cout<<Lucas(n+m,n)<<'\n';
    }
    return 0;
}

扩展卢卡斯定理

求 \(\binom{n}{m}\mod p\)，其中 \(p\) 不一定为质数。

设 \(p=p_1^{c_1}p_2^{c_2}...p_k^{c_k}\)，求出

\(\binom{n}{m}\mod p_1^{c_1}\)

\(\binom{n}{m}\mod p_2^{c_2}\)

\(\binom{n}{m}\mod p_k^{c_k}\)

之后用中国剩余定理进行合并即可。

现在求 \(\binom{n}{m}\mod p^k=\frac{n!}{m!(n-m)!}\mod p^k\)，但是无法通过求分母的逆元来求解，因为不一定有解。

考虑 \(a\) 对 \(p\) 有逆元的充要条件为 \(gcd(a,p)=1\)，于是将式子转化成 \(\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}p^{x-y-z}\mod p^k\)，其中 \(x\) 为 \(n!\) 中包含的 因子 \(p\) 的个数，\(y,z\) 同理，求出 \(\frac{n!}{p^x}\) 即可求逆元。

考虑如何求 \(\frac{n!}{p^x}\mod p^k,n!=(p*2p*3p...)(1*2*3...)\)，左半部分是 \(n\) 以内的 \(p\) 的倍数，右边是剩下的数。

\(1\) 到 \(n\) 中一共有 \(\left \lfloor \frac{n}{p} \right \rfloor\) 个 \(p\) 的倍数。

所以，\(n!=p^{\left \lfloor \frac{n}{p} \right \rfloor}*\left \lfloor \frac{n}{p} \right \rfloor!*(1*2*3...)\)，左边是 \(p\) 的个数，中间是所有 \(p\) 的倍数的系数，即这些倍数除以 \(p\) 的商，右边是剩余的数。

考虑对于右边剩下的部分 \(\prod_{i=1,i\not \equiv 0(\mod p)}^{n}i\)，而它对于 \(\mod p^k\) 是有循环节的。

即原式等于 \((\prod_{i=1,i\not \equiv 0(\mod p)}^{p^k}i)^{\left \lfloor \frac{n}{p^k} \right \rfloor}(\prod_{i=p^k\left \lfloor \frac{n}{p^k} \right \rfloor,i\not \equiv 0(\mod p)}^{n}i)\)，左边是 \(1\) 到 \(p^k\) 中所有无 \(p\) 因子的数的乘积，右边是余项的乘积。

回到 \(n!\) 的式子，前面的 \(p^{\left \lfloor \frac{n}{p} \right \rfloor}\) 是需要除掉的，但是 \(\left \lfloor \frac{n}{p} \right \rfloor!\) 中可能也有 \(p\) 因子，于是考虑递归。

定义 \(f(n)=\frac{n!}{p^x}\)，则 \(f(n)=f({\left \lfloor \frac{n}{p} \right \rfloor})(\prod_{i=1,i\not \equiv 0(\mod p)}^{p^k}i)^{\left \lfloor \frac{n}{p^k} \right \rfloor}(\prod_{i=p^k\left \lfloor \frac{n}{p^k} \right \rfloor,i\not \equiv 0(\mod p)}^{n}i)\)，递归边界 \(f(0)=1\).

回到最开始的式子，\(\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}p^{x-y-z}\mod p^k=\frac{f(n)}{f(m)f(n-m)}p^{x-y-z}\mod p^k\)，此时的 \(f(m)\) 一定与 \(p^k\) 互质，可以使用 \(exgcd\) 进行求解。

现在考虑计算最后剩下的 \(p^{x-y-z}\)，回到式子 \(n!=p^{\left \lfloor \frac{n}{p} \right \rfloor}*\left \lfloor \frac{n}{p} \right \rfloor!*(\prod_{i=1,i\not \equiv 0(\mod p)}^{p^k}i)^{\left \lfloor \frac{n}{p^k} \right \rfloor}(\prod_{i=p^k\left \lfloor \frac{n}{p^k} \right \rfloor,i\not \equiv 0(\mod p)}^{n}i)\)，需要计算的就是 \(p^{\left \lfloor \frac{n}{p} \right \rfloor}\) 且乘上的 \(\left \lfloor \frac{n}{p} \right \rfloor!\) 中的 \(p\) 因子。

于是得到了递推式子 \(g(n)=\left \lfloor \frac{n}{p} \right \rfloor+g(\left \lfloor \frac{n}{p} \right \rfloor)\)，递归边界 \(g(n)=0(\forall n<p)\)

所以得到最后的结论，\(\binom{n}{m}\mod p^k=\frac{f(n)}{f(m)f(n-m)}p^{g(n)-g(m)-g(n-m)}\mod p^k\)，最后再用中国剩余定理进行合并。

inline int CRT(int x,int p,int mod){/*合并某项该项同余方程，其中p为所有模数的lcm，mod为当前的模数*/
    return x*(p/mod)%p*inv(p/mod,mod)%p;
}
inline int fac(int n,int p,int pk){
    if(!n)return 1;/*递归边界f(0)=1*/
    int ans=1;
    for(int i=2;i<=pk;i++)if(i%p)(ans*=i)%=pk;/*p^k以内的不含p因子的数的乘积*/
    ans=po(ans,n/pk,pk);
    for(int i=2;i<=n%pk;i++)if(i%p)(ans*=i)%=pk;/*余项部分*/
    return ans*fac(n/p,p,pk)%pk;/*递归求解*/
}
inline int C(int n,int m,int p,int pk){/*计算C(n,m) mod p^k，其中p^k=pk*/
    int k=0;/*计算p^(x-y-z)，其中x-y-z=k*/
    for(int i=n;i;i/=p)k+=i/p;
    for(int i=m;i;i/=p)k-=i/p;
    for(int i=n-m;i;i/=p)k-=i/p;
    return fac(n,p,pk)*inv(fac(m,p,pk),pk)%pk*inv(fac(n-m,p,pk),pk)%pk*po(p,k,pk)%pk;
}
inline int exLucas(int n,int m,int p){
    int ans=0,tmp=p,pk;
    for(int i=2;i*i<=p;i++){
        if(tmp%i)continue;
        pk=1;
        while(tmp%i==0)pk*=i,tmp/=i;/*对p进行质因数分解得到pi^ki,其中pi=i,p^k=pk*/
        (ans+=CRT(C(n,m,i,pk),p,pk))%=p;/*计算并通过CRT合并C(n,m) mod p^k*/
    }
    if(tmp>1)(ans+=CRT(C(n,m,tmp,tmp),p,tmp))%=p;
    return ans;
}