【learning】 扩展lucas定理

首先说下啥是lucas定理：

$inom n m equiv inom {n\%P} {m\%P} imes inom{n/P}{m/P} pmod P$

借助这个定理，求$inom n m$时，若$P$较小，且$n,m$非常大时，我们就可以用这个定理要降低复杂度。

但是这个定理有一些限制，比如说要求$p$是质数，遇到一些毒瘤出题人不太好应对。

当$P$不是质数时，这时就要用到一个叫做扩展lucas定理的东西。

令$P=prod p_i^{k_i}$。

我们发现，如果对于每一个$p_i^{k_i}$，我们都求出$inom n m \% p_i^{k_i}$的值，我们就可以用$CRT$将它们合并，以得到最终的$inom n m$。

下面考虑如何求$inom n m \% p_i^{k_i}$。

首先考虑组合数的一个性质：

$inom n m =dfrac{n!}{m!(n-m)!}$

那么问题就可以归化为求n!及其逆元的问题

我们发现，我们可以考虑求出$n!$，$m!$的逆元，$(n-m)!$的逆元，然后就可以求出组合数了。

直接求的话，会发现$m!$和$(n-m)!$可能求不出逆元，因为$m!$可能会与$p_i^{k_i}$不互质。

我们定义$G(n,p_i)$表示$n!$中素因子$p_i$的个数，定义$F(n,p_i,k_i)=dfrac{n!}{p_i^{G(n,p_i)}}$。

则有：

$inom n m =dfrac{n!}{m!(n-m)!} equiv dfrac{F(n,p_i,k_i)p_i^{G(n,p_i)}}{F(m,p_i,k_i)p_i^{G(m,p_i)} imes  F(n-m,p_i,k_i) p_i^{G(n-m,p_i)}} pmod {p_i^{k_i}}$

考虑如何求函数$F$，我们显然有一种$O(n)$的求法：

$F(n,p_i,k_i)equiv prodlimits_{x=1,(x,p_i)=1}^{n}x imes F(n/p_i^{k_i},p_i,k_i) pmod {p_i^{k_i}}$

但是它依然是$O(n)$的

通过简单观察可以知道，求解F的连乘过程中有关于$p_i^{k_i}$的循环节，我们可以求出循环节的积，然后通过快速幂求解出前面若干个循环节的积的幂，最后乘上末尾非循环节部分的数。

举个例子：当$n=19,p=3,k=2$时：

$F(19,3,2)equiv 1 imes 2 imes 4 imes 5 imes 7 imes 8 imes 10 imes 11 imes 13 imes 16 imes 17 imes 19 imes F(6,3,2) pmod{p_i^{k_i}}$

$equiv (1 imes 2 imes 4 imes 5 imes 7 imes 8)^2 imes 19 imes pmod{p_i^{k_i}}$

根据这一个性质，我们得到：

$F(n,p_i,k_i)equiv left(prodlimits_{x=1,(n,p_i)=1}^{p_i^{k_i}} ight)^{left lfloor n/{p_i^{k_i}} ight floor}     F(n/p_i^{k_i},p_i,k_i) pmod{p_i^{k_i}}$

当$n=0$时，$F(n)=1$。

考虑如何求函数$G$，我们同样地采用递归的方式来搞，当$n>0$时，有：

$G(n,p_i)=lfloor frac{n}{p_i} floor +G(lfloor frac{n}{p_i} floor,p_i)$

当$n=0$时，$G$显然为$0$。

至此，我们求出了$inom n m \% p_i^{k_i}$。

我们求出了若干组这样的方程后，用CRT合并，就得到了最终的答案。

这种做法的复杂度也是非常地玄学，它是：

$O(sum p^klog(log_p n-k)+plog P)$

#include<bits/stdc++.h>
#define M 20000005
#define L long long
#define INF (1LL<<60)
using namespace std;

L pow_mod(L x,L k,const L MOD){L ans=1;for(;k;k>>=1,x=x*x%MOD) if(k&1) ans=ans*x%MOD; return ans;}

void exgcd(L a,L b,L &x,L &y){
    if(!b) {x=1; y=0; return;}
    exgcd(b,a%b,y,x);
    y-=a/b*x;
}
L inv(L a,L MOD){ L res1,res2; exgcd(a,MOD,res1,res2); return (res1+MOD)%MOD;}

L getp(L n,L p){L ans=0; for(;n;n/=p) ans=ans+n/p; return ans;}
L fac(L n,L p,L k){
    if(!n) return 1;
    L all=pow_mod(p,k,INF),mul=1,ans=1;
    for(L i=1;i<all;i++) if(i%p) mul=1LL*mul*i%all;
    ans=pow_mod(mul,n/all,all);
    for(L i=n%all;i;i--) if(i%p) ans=1LL*ans*i%all;
    return 1LL*ans*fac(n/p,p,k)%all;
}    

L get(L n,L m,L MOD){
    L c1=0,m1=1,p,up;
    for(p=2,up=sqrt(MOD);p<=up;p++) if(MOD%p==0){
        loop:;
        L k=(p>up),all=(p>up?p:1); 
        while(MOD%p==0) k++,MOD/=p,all*=p;
        L facn=fac(n,p,k);
        L facm=fac(m,p,k);
        L facnm=fac(n-m,p,k);
        L psum=getp(n,p)-getp(m,p)-getp(n-m,p);
        L c2=1LL*facn*inv(facm,all)%all*inv(facnm,all)%all*pow_mod(p,psum,all)%all;
        L mm=m1*all;
        L x=(1LL*inv(m1,all)*(c2-c1)%mm*m1+c1)%mm;
        m1=mm; c1=(x+m1)%m1;
    }
    if(MOD>1){p=MOD; MOD=1; goto loop;}
    return c1;
}        

main(){
    L z,y,p; cin>>z>>y>>p;
    cout<<get(z,y,p)<<endl;
}