NTTライブラリを持ってなかった。
http://agc005.contest.atcoder.jp/tasks/agc005_f
問題
N頂点の木を成す無向グラフが与えられる。
頂点群の部分集合Sに対し、f(S)はSを含む最小の連結した部分グラフの頂点数とする。
整数Kに対し、Sの選び方は通りある。その全ての選び方におけるf(S)の総和を求めたい。
K=1~Nそれぞれに対し、各総和を924844033で割った値求めよ。
解法
まずは各辺について考える。
辺の両側の先にa個・(N-a)個の頂点がつながっているとする。
Sとして頂点をK個選ぶとき、連結部分グラフにこの辺が含まれるには、辺の両側の頂点が1個以上選ばれた場合である。
よって逆に辺の片側だけの頂点が選ばれるケースを除くと、だけ解に寄与する。
また、上記値はあくまで辺の数だけで頂点数ではないので、頂点数は(辺の数+1)であることからだけさらに加える。
これを高速化することを考える。
b[i] := 各辺について、上記a・N-aを求めたとき、a=iまたは(N-a)=iを満たすa・(N-a)の数
とすると、上記総和はとなる。
N,Kが小さければ、各上記値の総和を答えればよいが、上記式の総和の部分はCombinationを前処理でO(1)で計算可能にしていてもO(N)なので、各Kに処理していては間に合わない。
各Kに対し上記総和の部分を早く求める事を考える。
= b[i]*i!, d[j]=1/(N-j)!と置くと、上記式はとなる。
こう式変形すると、総和の中は畳み込み演算の形なのでFFTで高速に処理できることがわかる。
幸い、924844033=2^21*441+1なので、w=5^441とするとw^(2^21) mod 924844033 = 1のためNTTで整数のまま計算可能。
ll mo=924844033; int N; vector<int> E[202020]; int C[202020]; int A[202020]; ll fact[202020]; ll combi(ll N_, ll C_) { const int NUM_=400001; static ll fact[NUM_+1],factr[NUM_+1],inv[NUM_+1]; if (fact[0]==0) { inv[1]=fact[0]=factr[0]=1; for (int i=2;i<=NUM_;++i) inv[i] = inv[mo % i] * (mo - mo / i) % mo; for (int i=1;i<=NUM_;++i) fact[i]=fact[i-1]*i%mo, factr[i]=factr[i-1]*inv[i]%mo; } if(C_<0 || C_>N_) return 0; return factr[C_]*fact[N_]%mo*factr[N_-C_]%mo; } void dfs(int cur,int pre) { C[cur]=1; FORR(r,E[cur]) if(r!=pre) dfs(r,cur),C[cur]+=C[r]; } ll modpow(ll a, ll n = mo-2) { ll r=1; while(n) r=r*((n%2)?a:1)%mo,a=a*a%mo,n>>=1; return r; } vector<ll> fft(vector<ll> v, bool rev=false) { int n=v.size(),i,j,m; 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]); } for(int m=2; m<=n; m*=2) { ll wn=modpow(5,(mo-1)/m); if(rev) wn=modpow(wn); for(i=0;i<n;i+=m) { ll w=1; for(int j1=i,j2=i+m/2;j2<i+m;j1++,j2++) { ll t1=v[j1],t2=w*v[j2]%mo; v[j1]=t1+t2; v[j2]=t1+mo-t2; while(v[j1]>=mo) v[j1]-=mo; while(v[j2]>=mo) v[j2]-=mo; w=w*wn%mo; } } } if(rev) { ll rv = modpow(n); FOR(i,n) v[i]=v[i]*rv%mo; } return v; } vector<ll> MultPoly(vector<ll> P,vector<ll> Q,bool resize=false) { if(resize) { int maxind=0,pi=0,qi=0,i; int s=2; FOR(i,P.size()) if(norm(P[i])) pi=i; FOR(i,Q.size()) if(norm(Q[i])) qi=i; maxind=pi+qi+1; while(s*2<maxind) s*=2; P.resize(s*2);Q.resize(s*2); } P=fft(P), Q=fft(Q); for(int i=0;i<P.size();i++) P[i]=P[i]*Q[i]%mo; return fft(P,true); } void solve() { int i,j,k,l,r,x,y; string s; cin>>N; FOR(i,N-1) { cin>>x>>y; E[x-1].push_back(y-1); E[y-1].push_back(x-1); } dfs(0,-1); FOR(x,N) if(C[x]!=N) A[C[x]]++,A[N-C[x]]++; vector<ll> P(1<<19),Q(1<<19); fact[0]=1; FOR(i,N+1) fact[i+1]=fact[i]*(i+1)%mo; for(i=0;i<=N;i++) { P[i] = A[i]*fact[i]; Q[N-i] = modpow(fact[i]); } vector<ll> R=MultPoly(P,Q); for(i=1;i<=N;i++) { ll ret=N*combi(N,i)%mo; ret -= R[N+i-0]*modpow(fact[i])%mo; cout<<(ret%mo+mo)%mo<<endl; } }
まとめ
素数サイコロと合成数サイコロも924844033の剰余で良ければもっと簡単だったのに…。