【HDU 6036】Division Game （NTT+数学）

#include <bits/stdc++.h>
#define rep(i,l,r) for(int i=l,ed=r;i<ed;++i)
#define mem(a,b) memset(a,b,sizeof(a))
typedef long long ll;
using namespace std;
const ll P = (235 << 22) + 1;
const int N = 1 << 18;
namespace Comb{
	const int N = 200005;
	ll fac[N] = {1, 1}, inv[N] = {1, 1}, f[N] = {1, 1};
	void init(ll M){
		for(int i = 2; i < N; i++) {
			fac[i] = fac[i - 1] * i % M;
			f[i] = (M - M / i) * f[M % i] % M;
			inv[i] = inv[i - 1] * f[i] % M;
		}
	}
	ll C(ll a, ll b) {
		if(b > a || b < 0)return 0;
		return fac[a] * inv[b] % P * inv[a - b] % P;
	}
};

ll qpow(ll a,ll b){
	ll ans=1;
	for(a%=P;b;b>>=1,a=a*a%P)if(b&1)ans=ans*a%P;
	return ans;
}

namespace NTT {
	ll w[N];
	int G=3;
	void ntt(ll *p,int n){
		for(int i=0,j=0;i<n;++i){
			if(i>j)swap(p[i],p[j]);
			for(int l=n>>1;(j^=l)<l;l>>=1);
		}
		for(int i=2;i<=n;i<<=1)
		for(int j=0,m=i>>1;j<n;j+=i)
			rep(k,0,m){
				ll b=w[n/i*k]*p[j+m+k]%P;
				p[j+m+k]=(p[j+k]-b+P)%P;
				p[j+k]=(p[j+k]+b)%P;
			}
	}
	void get_root(){
		static int pr[1000],cnt;
		int n=P-1;
		rep(i,2,sqrt(n)+0.5)
		if(n%i==0){
			pr[cnt++]=i;
			while(n%i==0)n/=i;
		}
		if(n>1)pr[cnt++]=n;
		rep(i,1,P){
			if(qpow(i,P-1)==1){
				bool fl=true;
				rep(j,0,cnt){
					if(qpow(i,(P-1)/pr[j])==1){
						fl=false;
						break;
					}
				}
				if(fl){
					G=i;break;
				}
			}
		}
	}
	void conv(int n,ll *x,ll *y){
		ll g=qpow(G,(P-1)/n);
		w[0]=1;
		rep(i,1,n)w[i]=w[i-1]*g%P;
		ntt(x,n);ntt(y,n);
		rep(i,0,n)x[i]=x[i]*y[i]%P;
		reverse(x+1,x+n);
		ntt(x,n);
	}
};
int m,k,p[20],e[20],w,n;
ll g[N],f[N],x[N],y[N];
int Cas;
int main(){
	Comb::init(P);
	//NTT::get_root();
	while(~scanf("%d%d",&m,&k)){
		w=1;
		rep(i,0,m){
			scanf("%d%d",p+i,e+i);
			w+=e[i];
		}
		
		rep(i,0,w)g[i]=1;
		rep(i,0,m)rep(j,0,w)g[j]=g[j]*Comb::C(e[i]+j-1,j-1)%P;
		
		mem(x,0);mem(y,0);
		rep(i,0,w){
			x[i]=(i&1?P-Comb::inv[i]:Comb::inv[i]);
			y[i]=g[i]*Comb::inv[i]%P;
		}
		
		for(n=1;n<w;n<<=1){}
		n<<=1;
		NTT::conv(n,x,y);
		
		ll invN=qpow(n,P-2);
		mem(f,0);
		rep(i,0,w)f[i]=x[i]*invN%P*Comb::fac[i]%P;
		
		printf("Case #%d:",++Cas);
		rep(i,1,k+1){
			ll tot=0;
			rep(j,0,w){
				tot+=qpow(f[j+1],i-1)*qpow(f[j],k-i+1)%P;
				if(tot>=P)tot-=P;
			}
			printf(" %lld",tot);
		}
		puts("");
	}
	return 0;
}