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

引子

求

[C_n^m ext{mod} p ]

不保证 (p) 是质数。

正文

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

前置知识

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

算法内容

将 (p) 用唯一分解定理分解，即

[p=prod p_i^{c_i} ]

若求出了

[{nchoose m} ext{mod} p_i^{c_i} ]

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

[{nchoose m} ext{mod} p^k (pin ext{prime})=cfrac{n!}{m!(n-m)!} ext{mod} p^k ]

但是显然 (m!) 和 ((n-m)!) 在 ( ext{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} ext{mod} p^k ]

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

将 (n!) 展开

[n!=1cdot2cdot3cdot...cdot n=(pcdot2pcdot3pcdot...)(1cdot2cdot...) ]

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

易知，在 (1) ~ (n) 中，有 (igglfloorcfrac{n}{p}igg floor) 个 (p) 的倍数。因此继续变形原式

[p^{lfloorfrac{n}{p} floor}(1cdot2cdot3cdot...)(1cdot2cdot...) ]

[p^{lfloorfrac{n}{p} floor} igg(igglfloorcfrac{n}{p}igg floorigg)! prodlimits_{i=1,i otequiv0pmod p}^n i ]

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

[p^{lfloorfrac{n}{p} floor} igg(igglfloorcfrac{n}{p}igg floorigg)! igg(prodlimits_{i=1,i otequiv0pmod p}^{p^k} iigg)^{lfloorfrac{n}{p^k} floor} igg(prodlimits_{i=p^klfloorfrac{n}{p^k} floor,i otequiv0pmod p}^{p^k} iigg) ]

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

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

[f(n)=figg(igglfloorcfrac{n}{p}igg floorigg) igg(prodlimits_{i=1,i otequiv0pmod p}^{p^k} iigg)^{lfloorfrac{n}{p^k} floor} igg(prodlimits_{i=p^klfloorfrac{n}{p^k} floor,i otequiv0pmod p}^{p^k} iigg) ]

且

[f(0)=0 ]

这样就可以 (O(nlog p)) 递推了。

回到最开始的式子

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

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

现在我完全明白了！

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

[n!=p^{lfloorfrac{n}{p} floor} igg(igglfloorcfrac{n}{p}igg floorigg)! igg(prodlimits_{i=1,i otequiv0pmod p}^{p^k} iigg)^{lfloorfrac{n}{p^k} floor} igg(prodlimits_{i=p^klfloorfrac{n}{p^k} floor,i otequiv0pmod p}^{p^k} iigg) ]

我们要求的部分就是 (p^{lfloorfrac{n}{p} floor}) 且乘上 (igg(igglfloorcfrac{n}{p}igg floorigg)!) 中的 (p) 因子。

因此我们得到递推式

[g(n)=igglfloorcfrac{n}{p}igg floor+g igg(igglfloorcfrac{n}{p}igg floorigg) ]

且

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

所以最后的结论就是

[{nchoose m} ext{mod} p^k (pin ext{prime})=cfrac{f(n)}{f(m)f(n-m)}cdot p^{g(n)-g(m)-g(n-m)} ext{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
",exLucas(n,m,p));
    return 0;
}

例题：礼物

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

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

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

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