【CF923E】Perpetual Subtraction
【CF923E】Perpetual Subtraction
by AmanoKumiko
Description
有一个整数\(x\)，范围在\([0,n]\)之间
初值为\(i\)的概率为\(p_i\)
然后会有\(m\)次操作
每次\(x\)会等概率变成\([0,x]\)中的整数
对于每个\(i\in [0,n]\)，求出\(m\)次操作后\(x\)为\(i\)的概率
对\(998244353\)取模
Input
第一行两个数\(n,m\)
第二行\(n+1\)个数读入\(p\)
Output
一行一个数表示答案
Sample Input
9 350
3 31 314 3141 31415 314159 3141592 31415926 314159265 649178508
Sample Output
822986014 12998613 84959018 728107923 939229297 935516344 27254497 413831286 583600448 442738326
Data Constraint
\(1\le n\le 10^5,1\le m\le 10^{18}\)
Solution
需要对积分比较敏感呢
令\(f_{i,j}\)表示\(i\)次操作后\(x\)为\(j\)的概率
那么有
\[f_{i,j}=\sum_{k=j}^n\frac{f_{i-1,k}}{k+1}\\
\]
用生成函数观察性质
\[\begin{align}
F_{i}(x)&=\sum_{j=0}^nf_{i,j}x^j\\
&=\sum_{j=0}^nx^j\sum_{k=j}^n\frac{f_{i-1,k}}{k+1}\\
&=\sum_{k=0}^n\frac{f_{i-1,k}}{k+1}\sum_{j=0}^kx^j\\
&=\sum_{k=0}^n\frac{f_{i-1,k}}{k+1}\frac{x^{k+1}-1}{x-1}\\
&=\frac{1}{x-1}\sum_{k=0}^nf_{i-1,k}\int_{1}^{x}t^{k}dt\\
&=\frac{1}{x-1}\int_{1}^{x}F_{i-1}(t)dt\\
\end{align}
\]
这里从\(1\)开始积分不好处理
我们希望的是\(\int_0^xf(t)dt=\int f(x)dx\)
那么设\(G_{i}(x)=F_{i}(x+1)\)同时\(G_{i}(x)=\sum_{j=0}^ng_{i,j}x^j\)
\[\begin{align}
G_{i}(x)&=\frac{1}{x}\int_{0}^{x}F_{i-1}(t+1)d(t+1)\\
&=\sum_{k=0}^{n}\frac{g_{i-1,k}x^{k}}{k+1}\\
\end{align}
\]
所以\(g_{i,j}=\frac{g_{i-1,j}}{j+1}\)
即\(g_{m,j}=\frac{g_{0,j}}{(j+1)^m}\)
众所周知，多项式平移可以用二项式定理展开变为卷积
\[\begin{align}
G_{0}(x)&=F_{0}(x+1)\\
&=\sum_{i=0}^np_i(x+1)^i\\
&=\sum_{i=0}^np_i\sum_{j=0}^i\binom{i}{j}x^j\\
&=\sum_{j=0}^n\frac{x^j}{j!}\sum_{i=j}^n\frac{i!p_i}{(i-j)!}
\end{align}
\]
\[\begin{align}
F_{m}(x)&=G_{m}(x-1)\\
&=\sum_{i=0}^ng_{m,i}(x-1)^i\\
&=\sum_{i=0}^ng_{m,i}\sum_{j=0}^i\binom{i}{j}x^j(-1)^{i-j}\\
&=\sum_{j=0}^n\frac{x^j}{j!}\sum_{i=j}^n\frac{g_{m,i}i!(-1)^{i-j}}{(i-j)!}
\end{align}
\]
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 500010
#define mo 998244353
#define LL long long
vector<int>Is;
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(){
Is.push_back(1);
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;
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)};
}
poly I(){return (poly){Is};}
int operator[](int x)const{
if(x<0||x>=val.size())return 0;
return val[x];
}
void NTT(int x){
F(i,0,sz()-1)if(rev[i]<i)swap(val[rev[i]],val[i]);
for(int mid=1;mid<sz();mid<<=1){
for(int i=0,gn=1;i<sz();i+=(mid<<1),gn=1){
F(j,0,mid-1){
int le=val[i+j],ri=1ll*gn*val[i+j+mid]%mo;
val[i+j]=mod(le+ri);val[i+j+mid]=mod(le-ri+mo);
gn=1ll*gn*(x==1?G1[mid]:G2[mid])%mo;
}
}
}
if(x==-1){F(i,0,sz()-1)val[i]=1ll*val[i]*inv[sz()]%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;
f.clear();g.clear();
g.val.push_back(mi(val[0],mo-2));
for(int i=2;i<Len*2;i<<=1){
f.rsz(i<<1);
g.rsz(i<<1);
BRT(i<<1);
F(j,0,i-1)f.val[j]=(j<val.size()?val[j]:0);
f.DFT();g.DFT();
F(j,0,(i<<1)-1)g.val[j]=1ll*g[j]*mod(mo+2-1ll*f[j]*g[j]%mo)%mo;
g.IDFT();
g.rsz(i);
}
return g.modxn(Len);
}
poly Ln(int Len){
return (deriv()*inver(Len)).integ().modxn(Len);
}
poly Exp(int Len){
poly f;
f.clear();
f.val.push_back(1);
for(int i=2;i<Len*2;i<<=1)f=(f*(I()-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;
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;
}
};
int n,p[N];
LL m;
int main(){
init();
scanf("%d%lld",&n,&m);
F(i,0,n)scanf("%d",&p[i]);
poly f,g,G0,Fm;
F(i,0,n)f.ins(ifac[i]);
F(i,0,n)g.ins(1ll*p[n-i]*fac[n-i]%mo);
f*=g;G0.rsz(n+1);
F(i,0,n)G0.val[i]=1ll*f[n-i]*ifac[i]%mo*mi(mi(i+1,m%(mo-1)),mo-2)%mo;
f.clear();g.clear();
F(i,0,n)f.ins((i&1)?mo-ifac[i]:ifac[i]);
F(i,0,n)g.ins(1ll*G0[n-i]*fac[n-i]%mo);
f*=g;Fm.rsz(n+1);
F(i,0,n)Fm.val[i]=1ll*f[n-i]*ifac[i]%mo;
F(i,0,n)printf("%d ",Fm[i]);
return 0;
}
