【CF923E】Perpetual Subtraction【生成函数】【积分推式子】【NTT卷积】

题意:有一个整数x∈[0,n]x\in[0,n]x[0,n],取iii的概率为pip_ipi。执行mmm次操作,每次把xxx等概率变成[0,x][0,x][0,x]中的一个整数,求操作完后等于每个数的概率。模998244353998244353998244353

n≤105,m≤1018n\leq 10^5,m\leq10^{18}n105,m1018

显然有一个dp式子

fi,j=∑k=jnfi−1,kk+1f_{i,j}=\sum_{k=j}^n\frac{f_{i-1,k}}{k+1}fi,j=k=jnk+1fi1,k

考虑生成函数

fi(x)=∑j=0n∑k=jn[xk]fi−1(x)k+1xjf_i(x)=\sum_{j=0}^n\sum_{k=j}^n\frac{[x^k]f_{i-1}(x)}{k+1}x^jfi(x)=j=0nk=jnk+1[xk]fi1(x)xj

=∑k=0n[xk]fi−1(x)k+1∑j=0kxj=\sum_{k=0}^n\frac{[x^k]f_{i-1}(x)}{k+1}\sum_{j=0}^kx^j=k=0nk+1[xk]fi1(x)j=0kxj

等比数列求和

=∑k=0n[xk]fi−1(x)k+1xk+1−1x−1=\sum_{k=0}^n\frac{[x^k]f_{i-1}(x)}{k+1}\frac{x^{k+1}-1}{x-1}=k=0nk+1[xk]fi1(x)x1xk+11

似乎推不动了

但是有奥妙重重的一步:

=∑k=0n[xk]fi−1(x)x−1xk+1−1k+1=\sum_{k=0}^n\frac{[x^k]f_{i-1}(x)}{x-1}\frac{x^{k+1}-1}{k+1}=k=0nx1[xk]fi1(x)k+1xk+11

看起来没啥用,但这样把右边写成了积分的形式

=∑k=0n[xk]fi−1(x)x−1∫1xtkdt=\sum_{k=0}^n\frac{[x^k]f_{i-1}(x)}{x-1}\int_1^xt^kdt=k=0nx1[xk]fi1(x)1xtkdt

然后盯着这个式子:

∑k=0n[xk]fi−1(x)∫1xtkdt\sum_{k=0}^n[x^k]f_{i-1}(x)\int_1^xt^kdtk=0n[xk]fi1(x)1xtkdt

ttt可能看不习惯,换个更通用的写法:

∑k=0n[xk]f(x)∫LRxkdx\sum_{k=0}^n[x^k]f(x) \int_L^Rx^kdxk=0n[xk]f(x)LRxkdx

相当于把多项式拆成单项式,每个单项式对区间[L,R][L,R][L,R]求定积分,最后还贴心地帮你乘上了系数

那不就是对多项式求积分吗

∫LRf(x)dx\int_L^Rf(x)dxLRf(x)dx

代回原式:

fi(x)=1x−1∫1xfi−1(t)dtf_i(x)=\frac{1}{x-1}\int_1^xf_{i-1}(t)dtfi(x)=x111xfi1(t)dt

然而这个积分下界是111,不好处理

然后又是奥妙重重的一步:

gi(x)=fi(x+1)g_i(x)=f_i(x+1)gi(x)=fi(x+1)

gi(x)=1x∫1x+1fi−1(t)dtg_i(x)=\frac{1}{x}\int_1^{x+1}f_{i-1}(t)dtgi(x)=x11x+1fi1(t)dt

=1x∫0xfi−1(t+1)dt=\frac{1}{x}\int_0^{x}f_{i-1}(t+1)dt=x10xfi1(t+1)dt

=1x∫0xgi−1(t)dt=\frac{1}{x}\int_0^{x}g_{i-1}(t)dt=x10xgi1(t)dt

那就gi−1g_{i-1}gi1的不定积分在xxx处的取值,容易得出

gi,j=1j+1gi−1,jg_{i,j}=\frac{1}{j+1}g_{i-1,j}gi,j=j+11gi1,j

这样如果知道g0g_0g0就可以轻松算出gmg_mgm

fffggg直接二项式定理暴拆再卷积一下就可以互相转换,问题解决

复杂度O(nlog⁡n)O(n\log n)O(nlogn)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#define MAXN 262144
using namespace std;
const int MOD=998244353;
typedef long long ll;
inline int read()
{
	int ans=0;
	char c=getchar();
	while (!isdigit(c)) c=getchar();
	while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
	return ans;
}
inline int qpow(int a,int p)
{
	int ans=1;
	while (p)
	{
		if (p&1) ans=(ll)ans*a%MOD;
		a=(ll)a*a%MOD;p>>=1;
	}
	return ans;
}
#define inv(x) qpow(x,MOD-2)
inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}
inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;}
int rt[2][24];
int r[MAXN],l,lim;
inline void init(){lim=1<<l;for (int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));}
inline void NTT(int* a,int type)
{
	for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);
	for (int L=0;L<l;L++)
	{
		int mid=1<<L,len=mid<<1;
		ll Wn=rt[type][L+1];
		for (int s=0;s<lim;s+=len)
			for (int k=0,w=1;k<mid;k++,w=w*Wn%MOD)
			{
				int x=a[s+k],y=(ll)w*a[s+mid+k]%MOD;
				a[s+k]=add(x,y),a[s+mid+k]=dec(x,y);
			}
	}
	if (type)
	{
		int t=inv(lim);
		for (int i=0;i<lim;i++) a[i]=(ll)a[i]*t%MOD;
	}
}
int fac[MAXN],finv[MAXN];
int a[MAXN];
int f[MAXN],g[MAXN];
int main()
{
	rt[0][23]=qpow(3,119);rt[1][23]=inv(rt[0][23]);
	for (int i=22;i>=0;i--)
	{
		rt[0][i]=(ll)rt[0][i+1]*rt[0][i+1]%MOD;
		rt[1][i]=(ll)rt[1][i+1]*rt[1][i+1]%MOD;
	}
	int n=read();
	ll m;
	scanf("%lld",&m);
	m%=MOD-1;
	for (int i=0;i<=n;i++) a[i]=read();
	fac[0]=1;
	for (int i=1;i<=n;i++) fac[i]=(ll)fac[i-1]*i%MOD;
	finv[n]=inv(fac[n]);
	for (int i=n-1;i>=0;i--) finv[i]=(ll)finv[i+1]*(i+1)%MOD;
	for (int i=0;i<=n;i++) f[i]=(ll)a[i]*fac[i]%MOD;
	for (int i=0;i<=n;i++) g[i]=finv[n-i];
	for (l=0;(1<<l)<=(n<<1);++l);
	init();
	NTT(f,0);NTT(g,0);
	for (int i=0;i<lim;i++) f[i]=(ll)f[i]*g[i]%MOD;
	NTT(f,1);
	for (int i=0;i<=n;i++) a[i]=(ll)f[i+n]*finv[i]%MOD;
	for (int i=0;i<=n;i++) a[i]=(ll)a[i]*inv(qpow(i+1,m))%MOD;
	for (int i=0;i<=n;i++) f[i]=(ll)a[i]*fac[i]%MOD;
	for (int i=0;i<=n;i++) g[i]=(((n-i)&1)? MOD-finv[n-i]:finv[n-i]);
	for (int i=n+1;i<lim;i++) f[i]=g[i]=0;
	NTT(f,0);NTT(g,0);
	for (int i=0;i<lim;i++) f[i]=(ll)f[i]*g[i]%MOD;
	NTT(f,1);
	for (int i=0;i<=n;i++) a[i]=(ll)f[i+n]*finv[i]%MOD;
	for (int i=0;i<=n;i++) printf("%d ",a[i]);
	return 0;
}
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值