【擴展lucas定理】

洛谷模板題面：https://www.luogu.org/problemnew/show/P4720
留坑講解……
先貼AC代碼：

#include<bits/stdc++.h>
using namespace std;
const int N =1000005;
#define rep(i,a,b) for(register int i=(a);i<=(b);++i)

typedef long long ll;
ll m,n;
int mod;
ll fac[N],inv[N];
ll ksm(ll x,ll y,ll M){
    ll aa=1ll;
    for(x%=M;y;y>>=1,x=(x*x)%M)if(y&1)aa=(aa*x)%M;
    return aa;
}
int p[N],pk[N],cnt;
ll sum,fak[22][N];
ll exgcd(ll x,ll y,ll &a,ll &b){
    if(!y){a=1,b=0;return x;}
    ll d=exgcd(y,x%y,b,a);
    b-=x/y*a;
    return d;
}
inline ll Inv(ll x,ll y){
    ll inv,rua;
    exgcd(x,y,inv,rua);
    return (inv+y)%y;
}
ll Fac(ll x,int i){
    if(x==0||x==1)return 1;
    return Fac(x/p[i],i)*ksm(fak[i][pk[i]-1],x/pk[i],pk[i])%pk[i]*fak[i][x%pk[i]]%pk[i];
}
ll ex_Lucas(ll x,ll y,int i){
    if(x<y)return 0;
    ll num=0;
    for(ll j=x;j;j/=p[i])
    num+=j/p[i];
    for(ll j=y;j;j/=p[i])num-=j/p[i];
    for(ll j=x-y;j;j/=p[i])num-=j/p[i];
    return Fac(x,i)*Inv(Fac(y,i),pk[i])%pk[i]*Inv(Fac(x-y,i),pk[i])*ksm(p[i],num,pk[i])%pk[i];
}
ll ans;
int main(){
    scanf("%lld%lld%d",&n,&m,&mod);
    int x=mod;
    for(int i=2;i*i<=mod;++i){
        if(x%i==0){
            p[++cnt]=i;
            pk[cnt]=1;
            while(x%i==0)x/=i,pk[cnt]*=i;
            sum=1;fak[cnt][0]=1;
            rep(j,1,pk[cnt]-1){
                if(j%p[cnt])sum=sum*j%pk[cnt];
                fak[cnt][j]=sum;
            }
        }
    }
    if(x!=1){
        ++cnt,p[cnt]=pk[cnt]=x;
        sum=1;fak[cnt][0]=1;
        rep(j,1,pk[cnt]-1){
            if(j%p[cnt])sum=sum*j%pk[cnt];
            fak[cnt][j]=sum;
        }
    }
    ll tmp;
    rep(i,1,cnt){
        tmp=ex_Lucas(n,m,i);
        ans=(ans+tmp*(mod/pk[i])%mod*Inv(mod/pk[i],pk[i])%mod)%mod;
    }
    printf("%lld\n",ans);
    return 0;
}