【学习笔记】扩展卢卡斯定理 exLucas

引子

求

\[C_n^m\ \text{mod}\ p
\]

不保证 $p$ 是质数。

正文

对于传统的 Lucas 定理，必须要求 $p$ 是质数才行。若 $p$ 不一定是质数，则需要扩展 Lucas 定理

前置知识

扩展欧几里得和中国剩余定理。

算法内容

将 $p$ 用唯一分解定理分解，即

\[p=\prod p_i^{c_i}
\]

若求出了

\[{n\choose m}\ \text{mod}\ p_i^{c_i}
\]

就可以用中国剩余定理合并答案了。那么此时我们要求的就是

\[{n\choose m}\ \text{mod}\ p^k\ (p\in \text{prime})=\cfrac{n!}{m!(n-m)!}\ \text{mod}\ p^k
\]

但是显然 $m!$ 和 $(n-m)!$ 在 $\text{mod}\ p^k$ 的意义下不一定存在逆元。

对于两个数 $a$ 和 $p$，存在逆元的充要条件是 $\gcd(a,p)=1$。

所以我们可以强行让他们互质。令 $x,y,z$ 分别为 $n!,m!,(n-m)!$ 中 $p$ 因子的个数，我们可以将原式转化得到

\[\cfrac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac {(n-m)!}{p^z}}\cdot p^{x-y-z}\ \text{mod} \ p^k
\]

这样我们就可以求逆元了，那么问题就转化成了对于每一个数 $n$ ，求出 $\cfrac{n!}{p^x}\ \text{mod}\ p^k$ 即可。

将 $n!$ 展开

\[n!=1\cdot2\cdot3\cdot...\cdot n=(p\cdot2p\cdot3p\cdot...)(1\cdot2\cdot...)
\]

其中左边的括号内是 $p$ 的所有倍数。

易知，在 $1$ ~ $n$ 中，有 $\bigg\lfloor\cfrac{n}{p}\bigg\rfloor$ 个 $p$ 的倍数。因此继续变形原式

\[p^{\lfloor\frac{n}{p}\rfloor}(1\cdot2\cdot3\cdot...)(1\cdot2\cdot...)
\]

\[p^{\lfloor\frac{n}{p}\rfloor}\ \bigg(\bigg\lfloor\cfrac{n}{p}\bigg\rfloor\bigg)!\ \prod\limits_{i=1,i\not\equiv0\pmod p}^n i
\]

然后你就发现他可以套娃了，把第三项拆成能整除之前和余数两个部分

\[p^{\lfloor\frac{n}{p}\rfloor}\ \bigg(\bigg\lfloor\cfrac{n}{p}\bigg\rfloor\bigg)!\ \bigg(\prod\limits_{i=1,i\not\equiv0\pmod p}^{p^k} i\bigg)^{\lfloor\frac{n}{p^k}\rfloor}\ \bigg(\prod\limits_{i=p^k\lfloor\frac{n}{p^k}\rfloor,i\not\equiv0\pmod p}^{p^k} i\bigg)
\]

也就是说上面这个式子是有循环节的，其中第三项是循环节 $1$ ~ $p$ 中非 $p$ 倍数的乘积，而第四项则是 $p$ 以后的。

前面这个 $p^{\lfloor\frac{n}{p}\rfloor}$ 是要除掉的，但是别忘了 $ \bigg(\bigg\lfloor\cfrac{n}{p}\bigg\rfloor\bigg)!$ 中可能也有 $p$ 因子。因此我们定义 $f(n)=\cfrac{n}{p^x}$，得到

\[f(n)=f\bigg(\bigg\lfloor\cfrac{n}{p}\bigg\rfloor\bigg)\ \bigg(\prod\limits_{i=1,i\not\equiv0\pmod p}^{p^k} i\bigg)^{\lfloor\frac{n}{p^k}\rfloor}\ \bigg(\prod\limits_{i=p^k\lfloor\frac{n}{p^k}\rfloor,i\not\equiv0\pmod p}^{p^k} i\bigg)
\]

且

\[f(0)=0
\]

这样就可以 $O(n\log p)$ 递推了。

回到最开始的式子

\[\cfrac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac {(n-m)!}{p^z}}\cdot p^{x-y-z}\ \text{mod} \ p^k=\cfrac{f(n)}{f(m)f(n-m)}\cdot p^{x-y-z}\ \text{mod}\ p^k
\]

此时 $f(m)$ 一定和 $p^k$ 互质，我们可以用扩展欧几里得求逆元来得到答案了。

但是还剩一个 $p^{x-y-z}$ 没算。比如，我们要算出现在 $f(n)$ 的 $x$。设 $g(n)=x$，再看这个式子

\[n!=p^{\lfloor\frac{n}{p}\rfloor}\ \bigg(\bigg\lfloor\cfrac{n}{p}\bigg\rfloor\bigg)!\ \bigg(\prod\limits_{i=1,i\not\equiv0\pmod p}^{p^k} i\bigg)^{\lfloor\frac{n}{p^k}\rfloor}\ \bigg(\prod\limits_{i=p^k\lfloor\frac{n}{p^k}\rfloor,i\not\equiv0\pmod p}^{p^k} i\bigg)
\]

我们要求的部分就是 $p^{\lfloor\frac{n}{p}\rfloor}$ 且乘上 $\bigg(\bigg\lfloor\cfrac{n}{p}\bigg\rfloor\bigg)!$ 中的 $p$ 因子。

因此我们得到递推式

\[g(n)=\bigg\lfloor\cfrac{n}{p}\bigg\rfloor+g \bigg(\bigg\lfloor\cfrac{n}{p}\bigg\rfloor\bigg)
\]

且

\[g(n)=0(\forall n<p)
\]

所以最后的结论就是

\[{n\choose m}\ \text{mod}\ p^k\ (p\in \text{prime})=\cfrac{f(n)}{f(m)f(n-m)}\cdot p^{g(n)-g(m)-g(n-m)}\ \text{mod}\ p^k
\]

然后用中国剩余定理合并答案就可以得到结果了！

代码

#include <bits/stdc++.h>

#define int long long

typedef long long ll;

const int maxn=1e3+10;

using namespace std;

inline int read(){

    int x=0;bool fopt=1;char ch=getchar();

    for(;!isdigit(ch);ch=getchar())if(ch=='-')fopt=0;

    for(;isdigit(ch);ch=getchar())x=(x<<3)+(x<<1)+ch-48;

    return fopt?x:-x;

}

inline int qpow(int x,int b,int p){

    int ans=1,base=x;

    while(b){

        if(b&1)ans=ans*base%p;

        base=base*base%p;

        b>>=1;

    }

    return ans;

}

void exgcd(int a,int b,int &x,int &y){

    if(!b)return x=1,y=0,void();

    exgcd(b,a%b,x,y);

    int z=x;x=y;y=z-a/b*y;

}

inline int inv(int a,int p){

    int x=0,y=0;

    exgcd(a,p,x,y);

    return (x+p)%p;

}

int F(int n,int p,int P){

    if(!n)return 1;

    int rou=1,rem=1;//rou表示循环节中的，rem表示余数部分的

    for(int i=1;i<=P;i++)

        if(i%p)rou=rou*i%P;

    rou=qpow(rou,n/P,P);

    for(int i=P*(n/P);i<=n;i++)

        if(i%p)rem=rem*(i%P)%P;

    return F(n/p,p,P)*rou%P*rem%P;

}

int G(int n,int p){

    if(n<p)return 0;

    return (n/p)+G(n/p,p);

}

inline int C(int n,int m,int p,int P){

    int x=F(n,p,P),y=inv(F(m,p,P),P),z=inv(F(n-m,p,P),P);

    int Pow=qpow(p,G(n,p)-G(m,p)-G(n-m,p),P);

    return x*y%P*z%P*Pow%P;

}

int tot;

int A[maxn],B[maxn];

inline void Init(int n,int m,int x){

    int M=sqrt(x);

    for(int i=2;i<=M;i++){

        if(x%i==0){

            int P=1;

            while(x%i==0){

                P*=i;x/=i;

            }

            A[++tot]=P;B[tot]=C(n,m,i,P);//算出每一个p^k意义下的答案

        }

    }

    if(x!=1){

        A[++tot]=x;B[tot]=C(n,m,x,x);

    }

}

inline int exLucas(int n,int m,int p){

    Init(n,m,p);

    int ans=0;

    for(int i=1;i<=tot;i++){

        int u=p/A[i],v=inv(u,A[i]);//普通的CRT合并答案

        ans=(ans+B[i]*u%p*v%p)%p;

    }

    return ans;

}

signed main(){

    int n=read(),m=read(),p=read();

    printf("%lld\n",exLucas(n,m,p));

    return 0;

}

例题：礼物

怎么这么多屑题叫礼物啊（恼）

显然答案是 $\prod\limits_{i=1}^mC_{n-\sum_{j=1}^{i-1}w_j}^{w_i}\pmod P$，然后用扩展卢卡斯求解就可以了。至于无解情况，只要看剩余的礼物数是否为负即可。