Codeforces 923E - Perpetual Subtraction（微积分+生成函数+推式子+二项式反演+NTT）

Codeforces 题目传送门 & 洛谷题目传送门

神仙题 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

首先考虑最朴素的 \(dp\)，设 \(dp_{z,i}\) 表示经过 \(z\) 次操作后剩下的数为 \(i\) 的概率，那么显然有 \(dp\) 转移方程 \(dp_{z,i}=\sum\limits_{j\ge i}dp_{z-1,j}·\dfrac{1}{j+1}\)。

边界条件 \(dp_{0,i}=p_i\)

直接递推显然不行，考虑优化，我们记 \(F_z(x)=\sum\limits_{i=0}^ndp_{z,i}x^i\)，我们不妨来探究一下 \(F_z(x)\) 有哪些性质：

\[\begin{aligned} F_z(x)&=\sum\limits_{i=0}^ndp_{z,i}x^i\\ &=\sum\limits_{i=0}^n\sum\limits_{j\ge i}dp_{z-1,j}\dfrac{1}{j+1}x^i\\ &=\sum\limits_{j=0}^ndp_{z-1,j}\dfrac{1}{j+1}\sum\limits_{i=0}^jx^i\\ &=\sum\limits_{j=0}^ndp_{z-1,j}\dfrac{1}{j+1}\dfrac{x^{j+1}-1}{x-1}\\ &=\dfrac{1}{x-1}\sum\limits_{j=0}^ndp_{z-1,j}\dfrac{x^{j+1}-1}{j+1} \end{aligned} \]

注意到

\[\dfrac{x^{j+1}}{j+1}=\int_0^xk^j\,\mathrm dk,\dfrac{1}{j+1}=\int_0^1k^j\,\mathrm dk \]

因此

\[\dfrac{x^{j+1}-1}{j+1}=\int_1^xk^j\,\mathrm dk \]

代入得

\[\begin{aligned} &\dfrac{1}{x-1}\sum\limits_{j=0}^ndp_{z-1,j}\dfrac{x^{j+1}-1}{j+1}\\ =&\dfrac{1}{x-1}\sum\limits_{j=0}^ndp_{z-1,j}\int_1^xk^j\,\mathrm dk\\ =&\dfrac{1}{x-1}\int_1^x\sum\limits_{j=0}^ndp_{z-1,j}k^j\,\mathrm dk\\ =&\dfrac{1}{x-1}\int_1^xF_{z-1}(k)\,\mathrm dk \end{aligned} \]

注意到这边积分下界是 \(1\)，并且左边还多出了个恶心的 \(\dfrac{1}{x-1}\)，非常棘手，因此考虑将这个 \(\dfrac{1}{x-1}\) 去掉，于是考虑设 \(G_z(x)=F_z(x+1)\)，那么：

\[\begin{aligned} G_z(x)&=F_z(x+1)\\ &=\dfrac{1}{x}\int_1^{x+1}F_{z-1}(k)\,\mathrm dk\\ &=\dfrac{1}{x}\int_0^{x}F_{z-1}(k+1)\,\mathrm dk\\ &=\dfrac{1}{x}\int_0^{x}G_{z-1}(k)\,\mathrm dk\\ \end{aligned} \]

噫！好！这下下界变成 \(0\) 了，左边也没有讨厌的 \(\dfrac{1}{x-1}\) 了！

接下来考虑设 \(G_z(x)=\sum\limits_{i=0}^ng_{z,i}x^i\)，那么根据 \(G_z(x)=\dfrac{1}{x}\int_0^xG_{z-1}(k)\,\mathrm dk\) 有

\[G_z(x)=\dfrac{1}{x}\int_0^xG_{z-1}(k)\,\mathrm dk=\dfrac{1}{x}\sum\limits_{i=1}^{n+1}\dfrac{g_{z-1,i-1}}{i}x^i=\sum\limits_{i=0}^{n}\dfrac{g_{z-1,i}}{i+1}x^i \]

因此 \(g_{z,i}=\dfrac{g_{z-1,i}}{i+1}\)。

这玩意儿显然可以快速幂加速，即 \(g_{m,i}=\dfrac{g_{0,i}}{(i+1)^m}\)。

也就是说接下来我们只需实现 \(F_0(x)\to G_0(x),G_m(x)\to F_m(x)\) 即可。

注意到

\[\begin{aligned} F_0(x+1)&=\sum\limits_{i=0}^np_i(x+1)^i\\ &=\sum\limits_{i=0}^np_i\sum\limits_{j=0}^i\dbinom{i}{j}x^j\\ &=\sum\limits_{j=0}^i(\sum\limits_{i=j}^n\dbinom{i}{j}p_i)x^j \end{aligned} \]

故

\[g_{0,i}=\sum\limits_{j=i}^n\dbinom{j}{i}p_j \]

拆组合数可得 \(g_{0,i}=\sum\limits_{j=i}^n\dfrac{j!}{i!(j-i)!}p_j=\dfrac{1}{i!}\sum\limits_{j=i}^nj!p_j·\dfrac{1}{(j-i)!}\)，按照 NTT 求差卷积的套路来即可。

\(G_m(x)\to F_m(x)\) 也同理，类似地有

\[g_{m,i}=\sum\limits_{j=i}^n\dbinom{j}{i}f_{m,j} \]

二项式反演一下可得

\[f_{m,i}=\sum\limits_{j=i}^n\dbinom{j}{i}(-1)^{j-i}g_{m,j} \]

和上面一样 NTT 求遍差卷积即可。

时间复杂度 \(n\log n\)。

const int MAXN=1e5;
const int MAXP=1<<18;
const int MOD=998244353;
const int pr=3;
const int ipr=(MOD+1)/3;
int n,f[MAXN+5],g[MAXN+5],h[MAXN+5];ll m;
int fac[MAXN+5],ifac[MAXN+5];
void init_fac(int n){
	for(int i=(fac[0]=ifac[0]=ifac[1]=1)+1;i<=n;i++) ifac[i]=1ll*ifac[MOD%i]*(MOD-MOD/i)%MOD;
	for(int i=1;i<=n;i++) fac[i]=1ll*fac[i-1]*i%MOD,ifac[i]=1ll*ifac[i-1]*ifac[i]%MOD;
}
int qpow(int x,int e){
	int ret=1;
	for(;e;e>>=1,x=1ll*x*x%MOD) if(e&1) ret=1ll*ret*x%MOD;
	return ret;
}
int rev[MAXP+5];
void NTT(vector<int> &a,int len,int type){
	int lg=31-__builtin_clz(len);
	for(int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
	for(int i=0;i<len;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=2;i<=len;i<<=1){
		int W=qpow((type<0)?ipr:pr,(MOD-1)/i);
		for(int j=0;j<len;j+=i){
			for(int k=0,w=1;k<(i>>1);k++,w=1ll*w*W%MOD){
				int X=a[j+k],Y=1ll*a[(i>>1)+j+k]*w%MOD;
				a[j+k]=(X+Y)%MOD;a[(i>>1)+j+k]=(X-Y+MOD)%MOD;
			}
		}
	}
	if(type==-1){
		int ivn=qpow(len,MOD-2);
		for(int i=0;i<len;i++) a[i]=1ll*a[i]*ivn%MOD;
	}
}
vector<int> conv(vector<int> a,vector<int> b){
	int LEN=1;while(LEN<a.size()+b.size()) LEN<<=1;
	a.resize(LEN,0);b.resize(LEN,0);NTT(a,LEN,1);NTT(b,LEN,1);
	for(int i=0;i<LEN;i++) a[i]=1ll*a[i]*b[i]%MOD;NTT(a,LEN,-1);
	return a;
}
int main(){
	scanf("%d%lld",&n,&m);m%=(MOD-1);init_fac(n);
	for(int i=0;i<=n;i++) scanf("%d",&f[i]);
	vector<int> A(n+1),B(n+1);
	for(int i=0;i<=n;i++) A[i]=ifac[i],B[n-i]=1ll*fac[i]*f[i]%MOD;
	vector<int> C=conv(A,B);
	for(int i=0;i<=n;i++){
		g[i]=1ll*ifac[i]*C[n-i]%MOD;
		g[i]=1ll*g[i]*qpow(qpow(i+1,MOD-2),m)%MOD;
	} iota(A.begin(),A.end(),0);iota(B.begin(),B.end(),0);
	for(int i=0;i<=n;i++){
		if(~i&1) A[i]=ifac[i];else A[i]=MOD-ifac[i];
		B[n-i]=1ll*fac[i]*g[i]%MOD;
	} C=conv(A,B);
	for(int i=0;i<=n;i++) printf("%d ",1ll*C[n-i]*ifac[i]%MOD);
	return 0;
}