如果这篇博客帮助到你，可以请我喝一杯咖啡~
CC BY-NC-SA 4.0 （除特别声明或转载文章外）
好题！教会了我相似对角化。
我们设 $f_{k,i}$ 表示 $k$ 次过后 $x=i$ 的概率，则有：
\[f_{k,i}=\sum\limits_{j=i}^n\frac{f_{k-1,j}}{j+1}\\\]令 $F_k$ 为列向量 \(\begin{bmatrix}f_{k,0}\\ \vdots\\f_{k,n}\end{bmatrix}\)，则有：
\[\begin{bmatrix}1&\frac{1}{2}&\cdots&\frac{1}{n+1}\\0&\frac{1}{2}&\cdots&\frac{1}{n+1}\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\frac{1}{n+1}\end{bmatrix}F_k=F_{k+1}\]设这个大方阵为 $D$，输入的列向量为 $p$。
我们要求的实际上就是 $D^mp$。
我们发现我们要求它的幂，考虑对角化这个矩阵。
简单地说，找到这个矩阵的所有特征向量，由于特征向量只会被拉伸，所以在特征基作为基向量的空间中，这个线性变换一定是对角矩阵，对角线上的值即为所有特征值，对角矩阵的幂是非常好计算的，复杂度为 $O(n\log n)$，再变换回一般的空间即可。设 $G$ 为将特征向量变为基向量的基变换矩阵，则有 $G^{-1}DG$ 是特征基下的该线性变换，所以我们有 $D^m=G(G^{-1}DG)^mG^{-1}$。
上三角矩阵的特征值就是对角线上的所有元素，所以有：
\[G^{-1}DG=\begin{bmatrix}1&0&\cdots&0\\0&\frac{1}{2}&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\frac{1}{n+1}\end{bmatrix}\]而特征向量是可以根据定义解得：
\[(D-\lambda I)x=0\]把 $\lambda=1,\cdots,\frac{1}{n+1}$ 代入，就可以解得 $G$ 为：
\[\begin{bmatrix}(-1)^{0}{0\choose 0}&(-1)^1{1\choose 0}&\cdots&(-1)^n{n\choose 0}\\0&(-1)^2{1\choose 1}&\cdots&(-1)^{n+1}{n\choose 1}\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&(-1)^{2n}{n\choose n}\end{bmatrix}\]一列就是一个特征向量。现在考虑求出 $G^{-1}$，打表发现：
\[G^{-1}=\begin{bmatrix}{0\choose 0}&{1\choose 0}&\cdots&{n\choose 0}\\0&{1\choose 1}&\cdots&{n\choose 1}\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&{n\choose n}\end{bmatrix}\]证明：
\[\sum\limits_{k=0}^n(-1)^{i+k}{k\choose i}{j\choose k}=[i=j]\] \[\begin{aligned} &\sum\limits_{k=0}^n(-1)^{i+k}{k\choose i}{j\choose k}\\ =&\sum\limits_{k=i}^j(-1)^{i+k}{j-i\choose k-i}{j\choose i}\\ =&{j\choose i}\sum\limits_{k=0}^{j-i}(-1)^{2i+k}{j-i\choose k}\\ =&{j\choose i}\sum\limits_{k=0}^{j-i}(-1)^{k}{j-i\choose k}\\ =&{j\choose i}[j-i=0]\\ =&[i=j] \end{aligned}\]现在问题变成给一个列向量左乘一个方阵，这个复杂度是 $O(n^2)$ 的，不能接受，但我们的方阵比较特殊。
先考虑 $Gp$ 的计算，设 $q=Gp$，则有
\[\begin{aligned} q_i=&\sum\limits_{k=i}^n(-1)^{i+k}{k\choose i}p_k\\ =&\sum\limits_{k=i}^n(-1)^{i+k}\frac{k!}{i!(k-i)!}p_k\\ =&\frac{(-1)^i}{i!}\sum\limits_{k=i}^n(-1)^{k}k!p_k\frac{1}{(k-i)!}\\ \end{aligned}\]容易发现这是一个减法卷积可以 $\text{NTT}$ 优化，左乘 $G^{-1}$ 的计算和这个区别不大，就不细说了，左乘 $(G^{-1}DG)^m$ 的计算就是给每一项乘上对角线上对应的那一个数即可。
那么这道题就结束了。
#include<cstdio>
#include<algorithm>
int const mod=998244353,g=3,gi=(mod+1)/g;
int p[400010],tmp[400010],r[400010],fac[400010],inv[400010],lim=1<<18,n,m;
long long nm;
int pow(int x,int y){
int res=1;
while(y){
if(y&1) res=1ll*res*x%mod;
x=1ll*x*x%mod;
y>>=1;
}
return res;
}
void ntt(int *a,int const &type){
for(int i=0;i<lim;i++)if(i<r[i])std::swap(a[i],a[r[i]]);
for(int mid=1;mid<lim;mid<<=1){
int rt=pow(type==1?g:gi,(mod-1)/(mid<<1));
for(int r=mid<<1,j=0;j<lim;j+=r)
for(int p=1,k=0;k<mid;k++,p=1ll*p*rt%mod){
int x=a[j+k],y=1ll*p*a[j+mid+k]%mod;
a[j+k]=(x+y)%mod,a[j+mid+k]=(x-y+mod)%mod;
}
}
if(type==-1)for(int p=pow(lim,mod-2),i=0;i<lim;i++)a[i]=1ll*a[i]*p%mod;
}
int main(){
for(int i=0;i<lim;i++)r[i]=(r[i>>1]>>1)|((i&1)*(lim>>1));
scanf("%d%lld",&n,&nm);
m=nm%(mod-1);
for(int i=0;i<=n;i++)scanf("%d",p+i);
fac[0]=1;
for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%mod;
inv[n]=pow(fac[n],mod-2);
for(int i=n;i;i--)inv[i-1]=1ll*inv[i]*i%mod;
for(int i=0;i<=n;i++)p[i]=1ll*p[i]*fac[i]%mod,tmp[i]=inv[i];
for(int i=0;i<=n;i++)if(i<n-i)std::swap(p[i],p[n-i]);
ntt(p,1),ntt(tmp,1);
for(int i=0;i<lim;i++)p[i]=1ll*p[i]*tmp[i]%mod;
ntt(p,-1);
for(int i=0;i<=n;i++)if(i<n-i)std::swap(p[i],p[n-i]);
for(int i=n+1;i<lim;i++)tmp[i]=p[i]=0;
for(int i=0;i<=n;i++)p[i]=1ll*p[i]*inv[i]%mod;
for(int i=0;i<=n;i++)p[i]=1ll*p[i]*pow(pow(i+1,mod-2),m)%mod;
for(int i=0;i<=n;i++)tmp[i]=inv[i],p[i]=1ll*p[i]*pow(mod-1,i)%mod*fac[i]%mod;
for(int i=0;i<=n;i++)if(i<n-i)std::swap(p[i],p[n-i]);
ntt(p,1),ntt(tmp,1);
for(int i=0;i<lim;i++)p[i]=1ll*p[i]*tmp[i]%mod;
ntt(p,-1);
for(int i=0;i<=n;i++)if(i<n-i)std::swap(p[i],p[n-i]);
for(int i=n+1;i<lim;i++)p[i]=0;
for(int i=0;i<=n;i++)p[i]=1ll*p[i]*pow(mod-1,i)%mod*inv[i]%mod;
for(int i=0;i<=n;i++)printf("%d ",p[i]);
return 0;
}