【模板】卢卡斯与拓展卢卡斯

LUCAS

$\left(\genfrac{}{}{0ex}{}{n}{m}\right)\equiv \left(\genfrac{}{}{0ex}{}{\lfloor \frac{n}{p}\rfloor }{\lfloor \frac{m}{p}\rfloor }\right)\cdot \left(\genfrac{}{}{0ex}{}{n\%p}{m\%p}\right)\mathrm{m}\mathrm{o}\mathrm{d}p$

---

EXLUCAS

设模数为$P(P=p_1^{k_1}{p}_2^{k_2}...p_n^{k^n}),p_i$为互不相同的质数。
首先分别求出$x_i=\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}$。
问题便转换成了求解：$x\equiv x_i\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}$。套用CRT解决。

那么问题在于如何快速求出所有的$\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}$。

当$k_i=1$时，套用LUCAS。

当$k_i>1$时，$\left(\genfrac{}{}{0ex}{}{n}{m}\right)=\frac{n!}{m!(n-m)!}$无法对质数的次幂快速取模，考虑分别对$n!,m!,(n-m)!$进行拆分取模，最终答案即为$n!\times inv(m!,p_i^{k_i})\times inv((n-m)!,p_i^{k_i})$($inv(x,p)$表示$x$在模$p$意义下的逆元，exgcd求解)。

拆分取模方式如下：

假设当前$p_i=3,k_i=2,P=3^2$，求解$19!\mathrm{m}\mathrm{o}\mathrm{d}P$。

$19!=(1\times 2\times 4\times 5\times 7\times ...\times 14\times 16\times 17\times 19)\times (3^6)\times (1\times 2\times 3\times 4\times 5\times 6)$

对于第一部分$1\times 2\times 4\times 5\times 7\times 8\equiv 10\times 11\times 13\times 14\times 16\times 17\mathrm{m}\mathrm{o}\mathrm{d}P$，等价于$(\prod _{i=1,i\nmid p_i}^{P-1}i{)}^{\lfloor \frac{n}{P}\rfloor }\mathrm{m}\mathrm{o}\mathrm{d}P$，剩余的数必然少于$P$个，而将右边的$\lfloor \frac{n}{p_i}\rfloor !$递归下去求解即可。

---

代码

#include<bits/stdc++.h>
#define debug(x) cerr<<#x<<":"<<x<<endl
using namespace std;
typedef long long ll;
const int N=5010;

ll n,m,p,q[N],a[N];
int tot;

namespace exlucas{
    
    inline int fp(ll x,ll y,ll p)
    {
        ll re=1;
        for(;y;y>>=1,x=x*x%p) 
         if(y&1) re=re*x%p;
        return re;
    }
    
    ll exgcd(ll n,ll m,ll &x,ll &y)
    {
        if(!m) 
          {x=1;y=0;return n;}
        ll re=exgcd(m,n%m,y,x);y-=(n/m)*x;
        return re;
    }
    
    inline ll inv(ll n,ll p)
    {
        ll x,y;
        exgcd(n,p,x,y);
        return (x%p+p)%p; 
    }
    
    inline ll C(ll n,ll m,ll p)
    {
        if(m>n) return 0;
        ll re=1,x,y,z,i;
        for(i=1;i<=m;++i){
            x=n-i+1;y=(int)inv(i,p);
            re=re*x%p*y%p; 
        }
        return re;
    }
    
    ll lucas(ll n,ll m,ll p)
    {
        if(m==0 || n==m) return 1;
        return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
    }
    
    ll cal(ll n,ll a,ll b,ll p)
    {
        if(!n) return 1;
        ll i,j,re=1;
        for(i=1;i<p;++i) if(i%a) re=re*i%p;
        re=fp(re,n/p,p);
        for(i=n/p*p+1;i<=n;++i) if(i%a) re=re*i%p;
        return re*cal(n/a,a,b,p)%p;
    }
    
    inline ll multilucas(ll n,ll m,ll a,ll b,ll p)
    {
        ll i,s=0;
        for(i=n;i;i/=a) s+=i/a;
        for(i=m;i;i/=a) s-=i/a;
        for(i=n-m;i;i/=a) s-=i/a;
        return fp(a,s,p)*cal(n,a,b,p)%p*inv(cal(m,a,b,p),p)%p*inv(cal(n-m,a,b,p),p)%p;
    }
    
    inline ll cal(ll n,ll m,ll p)
    {
    	ll i,j,b,c,d,ans,mod,x,y;
    	for(i=2;i*i<=p;++i)
    	 if(p%i==0){
    	 	q[++tot]=1;
    	 	for(j=0;p%i==0;p/=i) q[tot]*=i,j++;
            if(j==1) a[tot]=lucas(n,m,i);
            else a[tot]=multilucas(n,m,i,j,q[tot]);
    	 }
    	if(p>1){q[++tot]=p;a[tot]=lucas(n,m,p);}
        ans=a[1];mod=q[1];
        for(i=2;i<=tot;++i){
        	b=exgcd(mod,q[i],x,y);
        	c=a[i]-ans;
        	if(c%b) return -1LL;
        	d=q[i]/b;x=((c/b)%d*x%d+d)%d;
        	ans=ans+mod*x;
        	mod=mod/b*q[i];
        }
        return ans;
    }
}

int main(){
    scanf("%lld%lld%lld",&n,&m,&p);
    printf("%lld\n",exlucas::cal(n,m,p));
    return 0;
}