【JZOJ5740】【20190706】幻想世界

题目

小 $omega $ 想要进行烟火表演，她一开始有(n)颗彗星和(n)颗陨石

• 如果小 (omega) 有(i)颗彗星而没有陨石，那么她会消耗(i)颗彗星并得到(a_i)的火焰

• 如果小 (omega) 有(j)颗陨石而没有彗星，那么她会消耗(j)颗陨石并得到(a_j)的火焰

• 否则小 (omega) 有 (p) 的概率消耗一颗陨石，得到(alpha)的火焰，有 (q) 的概率消耗一颗彗星，得到(eta)的火焰

  有((1-p-q))的概率丢弃掉所有材料

  设(f_{i,j})表示小 (omega) 有 (i) 颗彗星， (j) 颗陨石期望得到的火焰数量

  求

  [sum_{i=0}^{n}sum_{j=0}^{n} h^{i(n+1)+j} f_{i,j} ]

$ n le 200000 $

题解

• 设(gamma = p alpha + q eta)，得到(f_{i,j})的转移式

  [egin{align} f_{i,j} = pf_{i-1,j} + qf_{i,j-1} + gamma end{align} ]

• 相当于向下贡献 $ imes p $ ，向右贡献 ( imes q) ，分别考虑(a_i,b_i,gamma)的贡献和

  [egin{align} ans &= \ &gamma sum_{i=1}^{n}sum_{j=1}^{n} sum_{x=i}^{n}sum_{y=j}^{n} (^{x-i+y-j}_{x-i})p^{x-i}q^{y-j} h^{x(n+1)+y}\ +&sum_{i=1}^{n}a_isum_{x=i}^{n}sum_{y=1}^{n}(^{x-i+y-1}_{x-i})p^{x-i}q^{y}h^{x(n+1)+y}\ +&sum_{j=1}^{n}b_jsum_{x=1}^{n}sum_{y=j}^{n}(^{x-1+y-j}_{y-j})p^{x}q^{y-j}h^{x(n+1)+y}\ +&sum_{i=1}^{n}a_ih^{i(n+1)} + b_{i}h^i \ end{align} ]

• 考虑前三项，交换求和顺序分别是：

  [egin{align} &sum_{x=0}^{n}sum_{y=0}^{n}(^{x+y}_x)p^xq^y(sum_{i=x+1}^{n}h^{i(n+1)})(sum_{j=y+1}^{n}h^{j}) \ &sum_{x=0}^{n}sum_{y=1}^{n}(^{x+y-1}_{x}) p^xq^y (sum_{i=x+1}^{n}a_{i-x}h^{i(n+1)})h^y \ &sum_{x=1}^{n}sum_{y=0}^{n}(^{x-1+y}_{y}) p^xq^y h^{x(n+1)} (sum_{j=y+1}^{n}b_{j-y}h^j) \ end{align} ]

• 拆开组合数可以直接看成关于(x,y)的卷积

• 第一个式子的系数直接预处理，第二个式子的系数可以预先做一遍卷积

• 话说我有些细节推错了好多遍．．．．．．．．

Code

//JZOJ锅掉了，暂时只知道这份代码能跑过大样例

#include<bits/stdc++.h>
#define ll long long 
#define mod 998244353
#define G 3
using namespace std;
const int N=1<<19;
int n,h,p,q,ta,tb,r,a[N],b[N],tx,ty,dx,dy,A[6][N],B[3][N],fac[N],inv[N],ny[N];
int hx[N],hy[N],sx[N],sy[N],sa[N],sb[N],rev[N],len,L,ans;
char gc(){
	static char*p1,*p2,s[1000000];
	if(p1==p2)p2=(p1=s)+fread(s,1,1000000,stdin);
	return(p1==p2)?EOF:*p1++;
}
int rd(){
	int x=0;char c=gc();
	while(c<'0'||c>'9')c=gc();
	while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+c-'0',c=gc();
	return x;
}
int pw(int x,int y){
	int re=1;if(y<0)y+=mod-1;
	for(;y;y>>=1,x=1ll*x*x%mod)
		if(y&1)re=1ll*re*x%mod;
	return re;
}
void inc(int&x,int y){x+=y;if(x>=mod)x-=mod;}
int pls(int x,int y){int re=x+y;return re>=mod?re-mod:re;}
int mis(int x,int y){int re=x-y;return re<0?re+mod:re;}
void ntt(int*F,int sg){
	for(int i=0;i<len;++i)if(i<rev[i])swap(F[i],F[rev[i]]);
	for(int i=1;i<len;i<<=1){
		int wn=pw(G,sg*(mod-1)/i/2);
		for(int j=0;j<len;j+=i<<1){
			int w=1;
			for(int k=0;k<i;++k,w=(ll)w*wn%mod){
				int x=F[j+k],y=(ll)w*F[j+k+i]%mod;
				F[j+k]=pls(x,y);F[j+k+i]=mis(x,y);
			}
		}
	}
	if(!~sg){
		int iv=pw(len,mod-2);
		for(int i=0;i<len;++i)F[i]=(ll)F[i]*iv%mod;
	}
}
int main(){
	freopen("dream.in","r",stdin);
	freopen("dream2.out","w",stdout);
	n=rd();h=rd();ta=rd();tb=rd();
	p=rd();p=(ll)p*pw(rd(),mod-2)%mod;
	q=rd();q=(ll)q*pw(rd(),mod-2)%mod;
	r=((ll)p*ta+(ll)q*tb)%mod;
	for(int i=1;i<=n;++i)a[i]=rd();
	for(int i=1;i<=n;++i)b[i]=rd();
	ty=pw(h,n+1),tx=pw(ty,n+1);
	dy=pw(h,mod-2),dx=pw(ty,mod-2);
	for(int i=n;i;--i){
		sx[i]=hx[i]=tx=(ll)tx*dx%mod;
		sy[i]=hy[i]=ty=(ll)ty*dy%mod;
		inc(sx[i],sx[i+1]);
		inc(sy[i],sy[i+1]);
		inc(ans,((ll)hx[i]*a[i]+(ll)hy[i]*b[i])%mod);
	}
	ny[1]=1;for(int i=2;i<=n<<1;++i)ny[i]=(ll)(mod-mod/i)*ny[mod%i]%mod;
	for(int i=fac[0]=inv[0]=1;i<=n<<1;++i){
		inv[i]=1ll*inv[i-1]*ny[i]%mod;
		fac[i]=1ll*fac[i-1]*i%mod;
	}

	for(L=0,len=1;len<=n<<1;len<<=1,L++);
	for(int i=0;i<len;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
		
	for(int i=1,tp=1,tq=1;i<=n;++i){
		A[2][i-1]=(ll)tp*sx[i]%mod*inv[i-1]%mod;
		A[3][i-1]=(ll)tq*sy[i]%mod*inv[i-1]%mod;
		tp=(ll)tp*p%mod;tq=(ll)tq*q%mod;
		A[0][i-1]=(ll)tp*hx[i]%mod*inv[i-1]%mod;
		A[1][i-1]=(ll)tq*hy[i]%mod*inv[i-1]%mod;
	}

	reverse(a,a+n+1);reverse(b,b+n+1);
	ntt(a,1);ntt(b,1);ntt(hx,1);ntt(hy,1);
	for(int i=0;i<len;++i){
		a[i]=(ll)a[i]*hx[i]%mod;
		b[i]=(ll)b[i]*hy[i]%mod;
	}
	ntt(a,-1);ntt(b,-1);//ntt(hx,-1);ntt(hy,-1);
	for(int i=0,tp=1,tq=1;i<=n;++i){
		A[4][i]=(ll)tp*a[n+i]%mod*inv[i]%mod;
		A[5][i]=(ll)tq*b[n+i]%mod*inv[i]%mod;
		tp=(ll)tp*p%mod;tq=(ll)tq*q%mod;
	}

	for(int i=0;i<6;++i)ntt(A[i],1);
	for(int i=0;i<len;++i){
	//	B[0][i]=(ll)A[2][i]*A[3][i]%mod;
	//	B[1][i]=(ll)A[4][i]*A[1][i]%mod;
	//	B[2][i]=(ll)A[0][i]*A[5][i]%mod;
		B[0][i]=((ll)A[2][i]*A[3][i]%mod*r+(ll)A[4][i]*A[1][i]+(ll)A[0][i]*A[5][i])%mod;
	}
	//for(int i=0;i<3;++i)ntt(B[i],-1);
	ntt(B[0],-1);
	for(int i=0;i<=n<<1;++i){
		//inc(ans , ((ll)r*B[0][i] + B[1][i] + B[2][i])%mod * fac[i] %mod);
		inc(ans , (ll)B[0][i]*fac[i]%mod);
	}
	cout<<ans<<endl;
	return 0;
}