[cf2015ICLFinalsDiv1J]Ceizenpok’s formula
阿新 • • 發佈:2017-10-02
pow print form final pac blog using cstring mes
題意:$C_n^m\% k$
解題關鍵:擴展lucas+中國剩余定理裸題
1 #include<algorithm> 2 #include<iostream> 3 #include<cstring> 4 #include<cstdio> 5 #include<cmath> 6 #include<cstdlib> 7 typedef long long ll; 8 using namespace std; 9 ll mod,n,m,x,y,module[10002],piset[10002],r[10002]; 10 11 ll mod_pow(ll x,ll n,ll p){ 12 ll res=1; 13 while(n){ 14 if(n&1) res=res*x%p; 15 x=x*x%p; 16 n>>=1; 17 } 18 return res; 19 } 20 21 ll extgcd(ll a,ll b,ll &x,ll &y){ 22 ll d=a; 23 if(b) d=extgcd(b,a%b,y,x),y-=a/b*x;24 else x=1,y=0; 25 return d; 26 } 27 28 ll inv(ll t,ll mod){ extgcd(t,mod,x,y);return (x+mod)%mod;} 29 30 ll multi(ll n,ll pi,ll pk){//求非互質的部分 31 if (!n) return 1; 32 ll ans=1; 33 for (ll i=2;i<=pk;i++) if(i%pi) ans=ans*i%pk; 34 ans=mod_pow(ans,n/pk,pk); 35 for(ll i=2;i<=n%pk;i++) if(i%pi) ans=ans*i%pk; 36 return ans*multi(n/pi,pi,pk)%pk; 37 } 38 39 40 ll exlucas(ll n,ll m,ll pi,ll pk){//組合數 c(n,m)mod pk=pi^k 41 if(m>n) return 0; 42 ll a=multi(n,pi,pk),b=multi(m,pi,pk),c=multi(n-m,pi,pk); 43 ll k=0; 44 for(ll i=n;i;i/=pi) k+=i/pi; 45 for(ll i=m;i;i/=pi) k-=i/pi; 46 for(ll i=n-m;i;i/=pi) k-=i/pi; 47 return a*inv(b,pk)%pk*inv(c,pk)%pk*mod_pow(pi,k,pk)%pk;//組合數求解完畢 48 } 49 50 ll crt(int n,ll *r,ll *m){ 51 ll M=1,ret=0; 52 for(int i=0;i<n;i++) M*=m[i]; 53 for(int i=0;i<n;i++){ 54 ll w=M/m[i]; 55 ret+=w*inv(w,m[i])*r[i]; 56 ret%=M; 57 } 58 return (ret+M)%M; 59 } 60 61 ll fz(ll n,ll *m,ll *piset){//分解質因子 62 ll num=0; 63 for (ll i=2;i*i<=n;i++){ 64 if(n%i==0){ 65 ll pk=1; 66 while(n%i==0) pk*=i,n/=i; 67 m[num]=pk; 68 piset[num]=i; 69 num++; 70 } 71 } 72 if(n>1) m[num]=n,piset[num]=n,num++; 73 return num; 74 } 75 76 ll excomb(ll n,ll m,ll p){ 77 ll num=fz(p,module,piset); 78 for(int i=0;i<num;i++){ 79 r[i]=exlucas(n,m,piset[i],module[i]); 80 } 81 return crt(num,r,module); 82 } 83 84 int main(){ 85 cin>>n>>m>>mod; 86 printf("%d",excomb(n,m,mod)); 87 return 0; 88 }
[cf2015ICLFinalsDiv1J]Ceizenpok’s formula