【学习笔记】扩展卢卡斯定理 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;
}

例题：礼物

一年一度的圣诞节快要来到了。每年的圣诞节小 E 都会收到许多礼物，当然他也会送出许多礼物。不同的人物在小 E 心目中的重要性不同，在小 E 心中分量越重的人，收到的礼物会越多。

小 E 从商店中购买了 \(n\) 件礼物，打算送给 \(m\) 个人，其中送给第 \(i\) 个人礼物数量为 \(w_i\)。请你帮忙计算出送礼物的方案数（两个方案被认为是不同的，当且仅当存在某个人在这两种方案中收到的礼物不同）。由于方案数可能会很大，你只需要输出模 \(P\) 后的结果。\(P\) 不保证是质数。若无法满足条件，输出 Impossible。

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

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