結果
問題 | No.215 素数サイコロと合成数サイコロ (3-Hard) |
ユーザー | kmjp |
提出日時 | 2015-05-30 02:10:41 |
言語 | C++11 (gcc 11.4.0) |
結果 |
WA
|
実行時間 | - |
コード長 | 7,060 bytes |
コンパイル時間 | 2,080 ms |
コンパイル使用メモリ | 179,024 KB |
実行使用メモリ | 12,468 KB |
最終ジャッジ日時 | 2024-07-06 12:31:54 |
合計ジャッジ時間 | 8,244 ms |
ジャッジサーバーID (参考情報) |
judge2 / judge3 |
(要ログイン)
ソースコード
#include <bits/stdc++.h> using namespace std; typedef signed long long ll; #undef _P #define _P(...) (void)printf(__VA_ARGS__) #define FOR(x,to) for(x=0;x<to;x++) #define FORR(x,arr) for(auto& x:arr) #define ITR(x,c) for(__typeof(c.begin()) x=c.begin();x!=c.end();x++) #define ALL(a) (a.begin()),(a.end()) #define ZERO(a) memset(a,0,sizeof(a)) #define MINUS(a) memset(a,0xff,sizeof(a)) //------------------------------------------------------- inline int mulmod(int a,int b,int mo) { int d,r; if(a==0 || b==0) return 0; __asm__("mull %4;" "divl %2" : "=d" (r), "=a" (d) : "g" (mo), "a" (a), "g" (b)); return r; } struct Kitamasa_fast { vector<int> P; static int ext_gcd(int p,int q,int& x, int& y) { // get px+qy=gcd(p,q) if(q==0) return x=1,y=0,p; int g=ext_gcd(q,p%q,y,x); y-=p/q*x; return g; } static int inv(int p,int q) { // return (1/p)%q ( p,q is co-prime) int xx,yy,g=ext_gcd(p,q,xx,yy); if(xx<0) xx+=q, yy-=p; return xx; } static int modpow(int a, int n, int mo) { int r=1; while(n) r=(n&1)?mulmod(r,a,mo):r,a=mulmod(a,a,mo),n>>=1; return r; } struct Mongo { ll Rmo; int S,R,R2,Ri,mo,Nd,mask; static int modpow(int a, int n, int mo) { int r=1; while(n) r=(n&1)?mulmod(r,a,mo):r,a=mulmod(a,a,mo),n>>=1; return r; } Mongo(int s,int m) { mo=m; S=s; R=1<<s; Rmo=1LL*R*mo; mask=R-1; R2=1LL*R*R%mo; Ri=modpow(R%mo,mo-2,mo); Nd = 0; int t=0, r=R, i=1; while(r>1) { if(t%2==0) t+=mo, Nd+=i; t/=2, r/=2, i*=2; } }; ll reduce (ll t) { t += ((t*Nd)&mask)*mo; t >>= S; if(t>=mo) t-=mo; return t; } ll conv(int t) { return reduce(1LL*t*R2);} ll rev(ll t) { return reduce(t);} ll mult(ll a,ll b) { return reduce(a*b);} // for convert int mult2(int a,int b) { return (int)reduce(reduce(1LL*a*b)*R2);} // for not converted ll modpow(ll a,ll n) { // a,r is converted ll r=conv(1); while(n) r=(n&1)?mult(r,a):r,a=mult(a,a),n>>=1; return r; } }; static vector<int> FMT(vector<int> v, int mo, bool rev=false) { // v.size()=2^k, mo = 2^22*k+1, int i,j,k,n=v.size(); // n=2^k; int pf=modpow(3,(mo-1)/n,mo); if(rev) pf=modpow(pf,mo-2,mo); for(i=0,j=1;j<n-1;j++) { for(int k=n>>1;k>(i^=k);k>>=1); if(i>j) swap(v[i],v[j]); } // mongo Mongo mongo(30,mo); vector<ll> V(n,0); FOR(i,n) V[i]=mongo.conv(v[i]); pf=mongo.conv(pf); for(int m=2; m<=n; m*=2) { ll ba=mongo.modpow(pf,n/m,mo); ll w=mongo.conv(1); FOR(i,m/2) { for(int j1=i,j2=i+m/2;j2<n;j1+=m,j2+=m) if(v[j2]) { ll t1=V[j1],t2=mongo.mult(w,V[j2]); V[j1]=t1+t2, V[j2]=t1-t2; if(V[j1]>=mongo.Rmo) V[j1]-=mongo.Rmo; if(V[j2]<0) V[j2]+=mongo.Rmo; } else { V[j2]=V[j1]; } w=mongo.mult(w,ba); } } if(rev) { ll rn=mongo.modpow(mongo.conv(n),mo-2); FOR(i,n) V[i]=mongo.mult(V[i],rn); } FOR(i,n) v[i]=mongo.rev(V[i]); /* for(int m=2; m<=n; m*=2) { int ba=modpow(pf,n/m,mo); int w=1; FOR(i,m/2) { for(int j1=i,j2=i+m/2;j2<n;j1+=m,j2+=m) if(v[j2]) { int t1=v[j1],t2=mulmod(w,v[j2],mo); v[j1]=t1+t2, v[j2]=t1-t2; if(v[j1]>=mo) v[j1]-=mo; if(v[j2]<0) v[j2]+=mo; } else { v[j2]=v[j1]; } w=mulmod(w,ba,mo); } } if(rev) { int rn=modpow(n,mo-2,mo); FOR(i,n) v[i]=mulmod(v[i],rn,mo); } */ return v; } static int CRT_garner(vector<pair<int,int> > V,int mo=1000000007) { int x,y,N=V.size(); static int invinv[3][3]; if(invinv[0][1]==0) { FOR(y,N) FOR(x,N) invinv[x][y]=inv(V[x].first,V[y].first); } FOR(y,N) FOR(x,y) { int k=V[y].second-V[x].second; if(k<0) k+=V[y].first; V[y].second = mulmod(k,invinv[x][y],V[y].first); } int ret=0; for(x=N-1;x>=0;x--) { ret = mulmod(ret,V[x].first,mo) + V[x].second; if(ret>=mo) ret -= mo; } return ret; } static vector<int> convol_sub(vector<int> a,vector<int> b, int mo,int t=99999999) { // mo = 2^k+1 int i,n=1; while(n<a.size()+b.size()-1) n*=2; a.resize(n); b.resize(n); a=FMT(a,mo); b=FMT(b,mo); FOR(i,n) a[i]=mulmod(a[i],b[i],mo); return FMT(a,mo,true); } static vector<int> convol(vector<int>& a,vector<int>& b, int mo,int t=99999999) { int mop[3]={0xA000001,0x1C000001,0x23800001}; auto x = convol_sub(a,b,mop[0],t); auto y = convol_sub(a,b,mop[1],t); auto z = convol_sub(a,b,mop[2],t); t=min(t,(int)x.size()); vector<int> R(t); vector<pair<int,int> > P{{mop[0],0},{mop[1],0},{mop[2],0}}; for(int i=0;i<t;i++) { // garner P[0].second = x[i]; P[1].second = y[i]; P[2].second = z[i]; R[i] = CRT_garner(P,mo); } return R; } static vector<int> mult(vector<int>& v,vector<int>& v2,vector<int>& D,vector<int>& id, int mo) { int k=v.size(),i; vector<int> res(k,0); vector<int> beta=convol(v,v2,mo); for(i=k-1;i<beta.size()&&i-(k-1)<k;i++) res[i-(k-1)]=beta[i]; beta.resize(k-1); vector<int> q=convol(beta,id,mo,k-1); q=convol(q,D,mo); for(i=k-1;i<q.size()&&i-(k-1)<k;i++) if((res[i-(k-1)]+=q[i])>=mo) res[i-(k-1)]-=mo; return res; } vector<int> getid(vector<int> D, int mo) { int t=1,k=D.size(); vector<int> id(1,1); while(t<=k) { t=min(2*t,k+1); vector<int> D2=D; D2.resize(t); vector<int> cur=convol(D2,id,mo,t); cur.resize(t); cur[0]+=2; if(cur[0]>=mo) cur[0]-=mo; id=convol(id,cur,mo,t); id.resize(t); } return id; } void calc(ll N, vector<int> D, int mo) { int n=D.size(); vector<int> p(n,0),v(n,0); p[0]=v[1]=1; reverse(ALL(D));reverse(ALL(p));reverse(ALL(v)); D.insert(D.begin(),mo-1); vector<int> id=getid(D,mo); while(N) { if(N%2) p=mult(p,v,D,id,mo); v=mult(v,v,D,id,mo); N/=2; } reverse(ALL(p)); P=p; } ll calc(ll N, vector<int> A, vector<int> D, int mo) { // A_K=A0*D0+A1*D1+A2*D2..+A_K-1*D_K-1 return A_N int i; ll ret=0; calc(N,D,mo); FOR(i,A.size()) ret += mulmod(A[i],P[i],mo); return (int)(ret%mo); } }; ll N; int P,C,M; int A[6]={2,3,5,7,11,13}; int B[6]={4,6,8,9,10,12}; int dp[2][320][4096]; ll dp2[8000]; ll mo=1000000007; void solve() { int i,j,k,l,r,x,y; string s; cin>>N>>P>>C; dp[0][0][0]=dp[1][0][0]=1; FOR(x,6) FOR(y,P) FOR(i,y*13+1) if((dp[0][y+1][i+A[x]] += dp[0][y][i])>=mo) dp[0][y+1][i+A[x]]-=mo; FOR(x,6) FOR(y,C) FOR(i,y*12+1) if((dp[1][y+1][i+B[x]] += dp[1][y][i])>=mo) dp[1][y+1][i+B[x]]-=mo; vector<int> p,q; FOR(i,4096) p.push_back(dp[0][P][i]), q.push_back(dp[1][C][i]); p=Kitamasa_fast::convol(p,q,mo); M=P*13+C*12; ll tot=0; if(N<=2*M) { dp2[0]=1; for(i=1;i<=2*M+2;i++) { FOR(x,M+1) if(i>=x) dp2[i] += dp2[i-x]*p[x]%mo; dp2[i] %= mo; } for(ll v=max(0LL,N-M);v<N;v++) { for(x=1;x<=M;x++) if(v+x>=N) tot += dp2[v]*p[x]%mo; tot %= mo; } } else { Kitamasa_fast kf; vector<int> A(M,1),V(M,0); FOR(i,M) V[i]=p[M-i]; tot = kf.calc(N+M-1,A,V,mo); } cout<<tot%mo<<endl; } int main(int argc,char** argv){ string s;int i; if(argc==1) ios::sync_with_stdio(false); FOR(i,argc-1) s+=argv[i+1],s+='\n'; FOR(i,s.size()) ungetc(s[s.size()-1-i],stdin); solve(); return 0; }