【JZOJ5740】【20190706】幻想世界

題目

小 $\omega $ 想要進行煙火表演，她一開始有\(n\)顆彗星和\(n\)顆隕石

• 若是小 \(\omega\) 有\(i\)顆彗星而沒有隕石，那麼她會消耗\(i\)顆彗星並獲得\(a_i\)的火焰

• 若是小 \(\omega\) 有\(j\)顆隕石而沒有彗星，那麼她會消耗\(j\)顆隕石並獲得\(a_j\)的火焰

• 不然小 \(\omega\) 有 \(p\) 的機率消耗一顆隕石，獲得\(\alpha\)的火焰，有 \(q\) 的機率消耗一顆彗星，獲得\(\beta\)的火焰

  有\((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 \beta\)，獲得\(f_{i,j}\)的轉移式
  \[ \begin{align} f_{i,j} \ = \ pf_{i-1,j} + qf_{i,j-1} + \gamma \end{align} \]

• 至關於向下貢獻 $\times p $ ，向右貢獻 \(\times q\) ，分別考慮\(a_i,b_i,\gamma\)的貢獻和
  \[ \begin{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_i\sum_{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_j\sum_{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} \]

• 考慮前三項，交換求和順序分別是：
  \[ \begin{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;
}