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 \times (n-m)! / p^c} * p^{a-b-c} \mod p^k$
然后就是$CRT$套上述这一堆东西
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<cstdlib> 5 #include<cmath> 6 #include<algorithm> 7 #include<queue> 8 #include<vector> 9 #include<map> 10 #include<set> 11 #define ll long long 12 #define db double 13 #define inf 2139062143 14 #define MAXN 200100 15 #define rep(i,s,t) for(register int i=(s),i##__end=(t);i<=i##__end;++i) 16 #define dwn(i,s,t) for(register int i=(s),i##__end=(t);i>=i##__end;--i) 17 #define ren for(register int i=fst[x];i;i=nxt[i]) 18 #define pb(i,x) vec[i].push_back(x) 19 #define pls(a,b) (a+b)%MOD 20 #define mns(a,b) (a-b+MOD)%MOD 21 #define mul(a,b) (1LL*(a)*(b))%MOD 22 using namespace std; 23 inline ll read() 24 { 25 ll x=0,f=1;char ch=getchar(); 26 while(!isdigit(ch)) {if(ch=='-') f=-1;ch=getchar();} 27 while(isdigit(ch)) {x=x*10+ch-'0';ch=getchar();} 28 return x*f; 29 } 30 ll T,n,m,MOD; 31 ll q_pow(ll a,ll t,ll p,ll res=1) 32 { 33 for(a%=p;t;t>>=1,(a*=a)%=p) 34 if(t&1) (res*=a)%=p;return res; 35 } 36 ll exgcd(ll a,ll b,ll &x,ll &y) 37 { 38 if(!b){x=1,y=0;return a;} 39 ll d=exgcd(b,a%b,y,x);y-=(a/b)*x;return d; 40 } 41 ll pw(ll n,ll p,ll pk) 42 { 43 if(!n) return 1;ll res=1; 44 rep(i,2,pk) if(i%p) (res*=i)%=pk; 45 res=q_pow(res,n/pk,pk); 46 rep(i,2,n%pk) if(i%p) (res*=i)%=pk; 47 return (res*pw(n/p,p,pk))%pk; 48 } 49 ll inv(ll n,ll p) {ll x,y;exgcd(n,p,x,y);return (x+p)%p;} 50 ll C(ll n,ll m,ll p,ll pk) 51 { 52 ll pn=pw(n,p,pk),pm=pw(m,p,pk),pz=pw(n-m,p,pk),sum=0; 53 for(ll i=n;i;i/=p) sum+=i/p;for(ll i=m;i;i/=p) sum-=i/p; 54 for(ll i=n-m;i;i/=p) sum-=i/p; 55 pm=inv(pm,pk),pz=inv(pz,pk); 56 return (((q_pow(p,sum,pk)*pn)%MOD*pm)%MOD*pz)%MOD; 57 } 58 void exlucas(ll n,ll m) 59 { 60 ll p=MOD,rs=p,k,ans=0,x,y; 61 rep(i,2,sqrt(MOD)) 62 { 63 k=1;while(rs%i==0) rs/=i,k*=i; 64 if(k!=1) (ans+=(inv(p/k,k)*p/k)%MOD*C(n,m,i,k)+MOD)%=MOD; 65 } 66 if(rs!=1) (ans+=(inv(p/rs,rs)*p/rs)%MOD*C(n,m,rs,rs)+MOD)%=MOD; 67 printf("%lld\n",ans); 68 } 69 int main() 70 { 71 n=read(),m=read(),MOD=read();exlucas(n,m); 72 }
