luogu 4720 【模板】扩展卢卡斯

题目大意：

求$C_n^m mod p$，p不一定为质数

思路：

首先可以将$p$分解为$p1^{a1}*p2^{a2}*...*pk^{ak}$，对于这些部分可以使用$CRT$合并

对于每个$p_i^{k_i}$，阶乘是存在循环的例如$19!$与模数$9$

$1*2*4*5*7*8$与$10*11*13*14*16*17$对答案的贡献一样，因此可以快速幂

对于剩下的部分因为很少可以暴力

对于求阶乘的部分 用这种方法求出循环节和剩余部分然后继续递归即可

求$C$的时候$C_n^m mod p^k= frac{n! / p^a}{m! / p^b imes (n-m)! / p^c} * p^{a-b-c} mod p^k$

然后就是$CRT$套上述这一堆东西

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<queue>
#include<vector>
#include<map>
#include<set>
#define ll long long
#define db double
#define inf 2139062143
#define MAXN 200100
#define rep(i,s,t) for(register int i=(s),i##__end=(t);i<=i##__end;++i)
#define dwn(i,s,t) for(register int i=(s),i##__end=(t);i>=i##__end;--i)
#define ren for(register int i=fst[x];i;i=nxt[i])
#define pb(i,x) vec[i].push_back(x)
#define pls(a,b) (a+b)%MOD
#define mns(a,b) (a-b+MOD)%MOD
#define mul(a,b) (1LL*(a)*(b))%MOD
using namespace std;
inline ll read()
{
    ll x=0,f=1;char ch=getchar();
    while(!isdigit(ch)) {if(ch=='-') f=-1;ch=getchar();}
    while(isdigit(ch)) {x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
ll T,n,m,MOD;
ll q_pow(ll a,ll t,ll p,ll res=1)
{
    for(a%=p;t;t>>=1,(a*=a)%=p)
        if(t&1) (res*=a)%=p;return res;
}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(!b){x=1,y=0;return a;}
    ll d=exgcd(b,a%b,y,x);y-=(a/b)*x;return d;
}
ll pw(ll n,ll p,ll pk)
{
    if(!n) return 1;ll res=1;
    rep(i,2,pk) if(i%p) (res*=i)%=pk;
    res=q_pow(res,n/pk,pk);
    rep(i,2,n%pk) if(i%p) (res*=i)%=pk;
    return (res*pw(n/p,p,pk))%pk;
}
ll inv(ll n,ll p) {ll x,y;exgcd(n,p,x,y);return (x+p)%p;}
ll C(ll n,ll m,ll p,ll pk)
{
    ll pn=pw(n,p,pk),pm=pw(m,p,pk),pz=pw(n-m,p,pk),sum=0;
    for(ll i=n;i;i/=p) sum+=i/p;for(ll i=m;i;i/=p) sum-=i/p;
    for(ll i=n-m;i;i/=p) sum-=i/p;
    pm=inv(pm,pk),pz=inv(pz,pk);
    return (((q_pow(p,sum,pk)*pn)%MOD*pm)%MOD*pz)%MOD;
}
void exlucas(ll n,ll m)
{
    ll p=MOD,rs=p,k,ans=0,x,y;
    rep(i,2,sqrt(MOD))
    {
        k=1;while(rs%i==0) rs/=i,k*=i;
        if(k!=1) (ans+=(inv(p/k,k)*p/k)%MOD*C(n,m,i,k)+MOD)%=MOD;
    }
    if(rs!=1) (ans+=(inv(p/rs,rs)*p/rs)%MOD*C(n,m,rs,rs)+MOD)%=MOD;
    printf("%lld
",ans);
}
int main()
{
    n=read(),m=read(),MOD=read();exlucas(n,m);
}