【LOJ6703】小Q的序列

【LOJ6703】小Q的序列

by AmanoKumiko

Description

给出一个长为\(n\)的序列\(c\)

求其所有非空子序列的权值和对\(998244353\)取模的结果

一个长为\(k\)序列\(a\)的权值为，\(\prod_{i=1}^k(a_i+i)\)

Input

第一行一个整数\(n\)

第二行\(n\)个数读入\(c\)

Output

一行一个数表示答案

Sample Input

3
1 2 3

Sample Output

90

Data Constraint

\(1\le n\le 10^5\)

Solution

直接上生成函数有点难以下手，先考虑\(dp\)

设\(f_{i,j}\)表示前\(i\)个数，选了\(j\)个的答案

那么有转移\(f_{i,j}=f_{i-1,j}+f_{i-1,j-1}(a_i+j)\)

但是会发现这里右边的\(j-1\)和\(a_i+j\)对上了

这使得我们上生成函数时难以利用导数合并（会得到一个这样鬼畜的柿子：\(F_{i}=(1+a_ix+x+x^2D)F_{i-1}(x)\)

不妨考虑令\(k=i-j\)

柿子变为\(f_{i,k}=f_{i-1,k-1}+f_{i-1,k}(a_i+i-k)\)

Sol1

\[\begin{align} F_i(x)&=xF_{i-1}(x)+\sum_{j=0}^{+\infty}f_{i-1,j}x^j(a_i+i-j)\\ &=xF_{i-1}(x)+(a_i+i)F_{i-1}(x)-x\sum_{j=1}^{+\infty}f_{i-1,j}jx^{j-1}\\ &=xF_{i-1}(x)+(a_i+i)F_{i-1}(x)-xF_{i-1}(x)'\\ \end{align} \]

观察：\((fe^{-x})'=f'e^{-x}-fe^{-x}\)

于是有

\[\begin{align} F_i(x)e^{-x}&=x(F_{i-1}(x)-F_{i-1}(x)')e^{-x}+(a_i+i)F_{i-1}(x)e^{-x}\\ F_i(x)e^{-x}&=(a_i+i)F_{i-1}(x)e^{-x}-x(F_{i-1}(x)e^{-x})'\\ \end{align} \]

令\(G_i(x)=F_{i}(x)e^{-x}\)

\[\begin{align} G_i(x)&=(a_i+i)G_{i-1}(x)-xG_{i-1}'\\ \sum g_{i,j}x^j&=(a_i+i)\sum g_{i-1,j}x^j-x\sum g_{i-1,j}x^{j-1}j\\ \sum g_{i,j}x^j&=(a_i+i)\sum g_{i-1,j}x^j-\sum g_{i-1,j}x^{j}j\\ g_{i,j}&=g_{i-1,j}(a_i+i-j)\\ g_{n,k}&=g_{0,k}\prod_{i=1}^n(a_i+i-k)\\ \end{align} \]

那么我们令\(P(x)=\prod_{i=1}^n(a_i+i-x)\)

分治\(NTT\)算出来最后多点求值就行了

Sol2

由于这题是对所有项系数求和，

于是有了一个不用多点求值的常数更小的做法

\[\begin{align} F_i(x)&=xF_{i-1}(x)+(a_i+i)F_{i-1}(x)-xF_{i-1}(x)'\\ &=(x+a_i+i-xD)F_{i-1}(x)\\ \end{align} \]

（\(D\)是微分算子\(\frac{d}{dx}\)）

那么我们还是设\(P(x)=\prod_{i=1}^n(a_i+i+x)\)，最后把\((x-xD)\)代进去即可

（这个地方的可行性可能对三观的冲击不小，不过还是应该好好理解一下。。。

那么我们要求的就是\(G_i(x)=(x-xD)^i\)的各项系数

考虑递推

\[\begin{align} G_i(x)&=(x-xD)^i\\ &=xG_{i-1}(x)-xDG_{i-1}(x)\\ \sum g_{i,j}x^j&=x\sum g_{i-1,j}x^j-\sum g_{i-1,j}x^{j}j\\ g_{i,j}&=-jg_{i-1,j}+g_{i-1,j-1} \end{align} \]

对比斯特林数\(S(n,m)=mS(n-1,m)+S(n-1,m-1)\)

我们可以写出\(g_{i,j}\)的组合意义：

将\(i\)个不同的元素放进\(j\)个集合中，每个集合若有\(x\)个元素，则权值为\((-1)^{x+1}\)

求所有方案下集合权值的乘积和

于是我们可以写出一个集合的\(EGF\)：\(\sum_{i=1}^{+\infty}\frac{(-1)^{i+1}x^i}{i!}\)

即\(-e^{-x}+1\)

那么我们的\((x-xD)^i\)求的是一行的系数和

同样地，仿照斯特林数，我们给\(-e^{-x}+1\)求个\(\exp\)即可

Code

#include<bits/stdc++.h>
using namespace std;
#define F(i,a,b) for(int i=a;i<=b;i++)
#define Fd(i,a,b) for(int i=a;i>=b;i--)
#define N 600010
#define mo 998244353
#define LL long long
#define ULL unsigned long long

int rev[N],G1[N],G2[N],fac[N],ifac[N],inv[N];

int mod(int x){return x>=mo?x-mo:x;}

int mi(int x,int y){
	if(!y)return 1;
	if(y==1)return x;
	return y%2?1ll*x*mi(1ll*x*x%mo,y/2)%mo:mi(1ll*x*x%mo,y/2);
}

void init(){
	fac[0]=ifac[0]=1;
	F(i,1,N-10)fac[i]=1ll*fac[i-1]*i%mo,inv[i]=(i==1?1:1ll*mo/i*mod(mo-1ll*inv[mo%i]%mo)%mo);
	ifac[N-10]=mi(fac[N-10],mo-2);
	Fd(i,N-11,1)ifac[i]=1ll*ifac[i+1]*(i+1)%mo;
	for(int l=1;l<=N-10;l<<=1)G1[l]=mi(3,(mo-1)/(l*2)),G2[l]=mi(G1[l],mo-2);
}

void BRT(int x){F(i,0,x-1)rev[i]=(rev[i>>1]>>1)|((i&1)?(x>>1):0);}

struct poly{
	vector<int>val;
	poly(int x=0){if(x)val.push_back(x);}
	poly(const vector<int>&x){val=x;}
	void Rev(){reverse(val.begin(),val.end());}
	void ins(int x){val.push_back(x);}
	void clear(){vector<int>().swap(val);}
	int sz(){return val.size();}
	void rsz(int x){val.resize(x);}
	void shrink(){for(;sz()&&!val.back();val.pop_back());}
	poly modxn(int x){
		if(val.size()<=x)return poly(val);
		else return poly(vector<int>(val.begin(),val.begin()+x));
	}
	int operator[](int x)const{
		if(x<0||x>=val.size())return 0;
		return val[x];
	}
	void NTT(int x){
		static ULL f[N],w[N];
		w[0]=1;
		F(i,0,sz()-1)f[i]=(((LL)mo<<5)+val[rev[i]])%mo;
		for(int mid=1;mid<sz();mid<<=1){
			int tmp=(x==1?G1[mid]:G2[mid]);
			F(i,1,mid-1)w[i]=w[i-1]*tmp%mo;
			for(int i=0;i<sz();i+=(mid<<1)){
				F(j,0,mid-1){
					int t=w[j]*f[i|j|mid]%mo;
					f[i|j|mid]=f[i|j]+mo-t;f[i|j]+=t;
				}
			}
			if(mid==(1<<10)){F(i,0,sz()-1)f[i]%=mo;};
		}
		if(x==-1){int tmp=inv[sz()];F(i,0,sz()-1)val[i]=f[i]%mo*tmp%mo;}
		else{F(i,0,sz()-1)val[i]=f[i]%mo;}
	}
	void DFT(){NTT(1);}
	void IDFT(){NTT(-1);}
	friend poly operator*(poly x,poly y){
		if(x.sz()<30||y.sz()<30){
			if(x.sz()>y.sz())swap(x,y);
			poly ret;
			ret.rsz(x.sz()+y.sz());
			F(i,0,ret.sz()-1){
				for(int j=0;j<=i&&j<x.sz();j++)
					ret.val[i]=mod(ret.val[i]+1ll*x[j]*y[i-j]%mo);
			}
	//		ret.shrink();
			return ret;
		}
		int l=1;
		while(l<x.sz()+y.sz()-1)l<<=1;
		x.rsz(l);y.rsz(l);BRT(l);
		x.DFT();y.DFT();
		F(i,0,l-1)x.val[i]=1ll*x[i]*y[i]%mo;
		x.IDFT();
	//	x.shrink();
		return x;
	}
	friend poly operator+(poly x,poly y){
		poly ret;
		ret.rsz(max(x.sz(),y.sz()));
		F(i,0,ret.sz()-1)ret.val[i]=mod(x[i]+y[i]);
		return ret;
	}
	friend poly operator-(poly x,poly y){
		poly ret;
		ret.rsz(max(x.sz(),y.sz()));
		F(i,0,ret.sz()-1)ret.val[i]=mod(x[i]-y[i]+mo);
		return ret;
	}
	poly &operator*=(poly x){return (*this)=(*this)*x;}
	poly &operator+=(poly x){return (*this)=(*this)+x;}
	poly &operator-=(poly x){return (*this)=(*this)-x;}
	poly deriv(){
		poly f;
		f.rsz(sz()-1);
		F(i,0,sz()-2)f.val[i]=1ll*(i+1)*val[i+1]%mo;
		return f;
	}
	poly integ(){
		poly f;
		f.rsz(sz()+1);
		F(i,1,sz())f.val[i]=1ll*val[i-1]*inv[i]%mo;
		return f;
	}
	poly inver(int Len){
		poly f,g,res(mi(val[0],mo-2));
		for(int i=1;i<Len;){
			i<<=1;f.rsz(i);g.rsz(i);BRT(i);
			F(j,0,i-1)f.val[j]=(*this)[j],g.val[j]=res[j];
			f.DFT();g.DFT();
			F(j,0,i-1)f.val[j]=1ll*f[j]*g[j]%mo;
			f.IDFT();
			F(j,0,(i>>1)-1)f.val[j]=0;
			f.DFT();
			F(j,0,i-1)f.val[j]=1ll*f[j]*g[j]%mo;
			f.IDFT();
			res.rsz(i);
			F(j,i>>1,i-1)res.val[j]=mod(mo-f[j]);
		}
		return res.modxn(Len);
	}
	poly Sqrt(int Len){
		poly f,g,res(1);
		for(int i=1;i<Len;){
			i<<=1;
			f=res;
			f.rsz(i>>1);BRT(i>>1);
			f.DFT();
			F(j,0,(i>>1)-1)f.val[j]=1ll*f[j]*f[j]%mo;
			f.IDFT();
			F(j,0,i-1)f.val[j%(i>>1)]=mod(f[j%(i>>1)]+mo-(*this)[j]);
			g=(2*res).inver(i>>1);f=(f*g).modxn(i>>1);f.rsz(i);
			F(j,i>>1,i-1)f.val[j]=f[j-(i>>1)];
			F(j,0,(i>>1)-1)f.val[j]=0;
			res-=f;
		}
		return res.modxn(Len);
	}
	poly Ln(int Len){
		return (deriv()*inver(Len)).integ().modxn(Len);
	}
	poly Exp(int Len){
		poly f(1);
		for(int i=2;i<Len*2;i<<=1)f=(f*(1-f.Ln(i)+modxn(i))).modxn(i);
		return f.modxn(Len);
	}
	poly Pow(int Len,int k){
		poly f;
		f.clear();
		int tail=0;
		while(val[tail]==0&&tail<sz())tail++;
		if(tail>=sz())return f;
		if(tail*k>=Len)return f;
		f.rsz(Len);
		int Mul=mi(val[tail],mo-2);
		F(i,0,min(Len-1,sz()-tail-1))f.val[i]=1ll*val[i+tail]*Mul%mo;
		Mul=mi(val[tail],k);
		f=f.Ln(Len);
		F(i,0,Len-1)f.val[i]=1ll*f[i]*(k%mo)%mo;
		f=f.Exp(Len);
		Fd(i,Len-1,tail*k)f.val[i]=1ll*f[i-tail*k]*Mul%mo;
		F(i,0,tail*k-1)f.val[i]=0;
		return f;
	}
};


poly f[N];

int n,c[N],ans;

poly solve(int l,int r){
	if(l==r)return f[l];
	int mid=l+r>>1;
	return solve(l,mid)*solve(mid+1,r);
}

int main(){
	init();
	scanf("%d",&n);
	F(i,1,n){
		scanf("%d",&c[i]);
		f[i].ins(c[i]+i);
		f[i].ins(1);
	}
	poly g=solve(1,n),h;
	F(i,0,n)h.ins((i&1)?ifac[i]:mo-ifac[i]);
	h+=1;
	h=h.Exp(n+1);
	F(i,0,n)ans=mod(ans+1ll*g[i]*h[i]%mo*fac[i]%mo);
	printf("%d",mod(ans-1+mo));
	return 0;
}