【LOJ6703】小Q的序列

阿新 • • 發佈：2022-04-08

【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

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;
}

【LOJ6703】小Q的序列