\(\mathcal{Description}\)
Link.
有一个整数 \(x\in[0,n]\)，初始时以 \(p_i\) 的概率取值 \(i\)。进行 \(m\) 轮变换，每次均匀随机取整数 \(r\in[0,x]\)，令 \(x\leftarrow r\)。求变换完成后 \(x=i~(i=0..n)\) 的概率。答案模 \(998244353\)。
\(\mathcal{Solution}\)
令向量 \(\boldsymbol p\) 为此时 \(x\) 的取值概率，显然一次变换是对 \(\boldsymbol p\) 的线性变换 \(A\)，有
\[A=\begin{bmatrix}
1&\frac{1}{2}&\cdots&\frac{1}{n}&\frac{1}{n+1}\\
&\frac{1}{2}&\cdots&\frac{1}{n}&\frac{1}{n+1}\\
&&\ddots&\vdots&\vdots\\
&&&\frac{1}{n}&\frac{1}{n+1}\\
&&&&\frac{1}{n+1}
\end{bmatrix}.
\]
我们的目标即为求出 \(A^m\boldsymbol p\)。注意到 \(A\) 的特征值很明显——\(\lambda_i=\frac{1}{i+1},i\in[0,n]\)，所以可以考虑将其对角化来加速矩阵幂的计算。手算一下 \(\lambda_i\) 所对应的特征向量 \(\boldsymbol v_i\)，发现一组特解
\[\boldsymbol v_i=\begin{bmatrix}
\binom{i}{0}\\
-\binom{i}{1}\\
\binom{i}{2}\\
\vdots\\
(-1)^n\binom{i}{n}
\end{bmatrix},
\]
那么 \(V=\begin{bmatrix}\boldsymbol v_0&\boldsymbol v_1&\cdots&\boldsymbol v_n\end{bmatrix}\) 有
\[V_{ij}=(-1)^i\binom{j}{i}.
\]
尝试对 \(V\) 求逆，由于
\[\begin{aligned}
V_{ij}^2 &= \sum_{k=0}^{n-1}(-1)^i\binom{k}{i}\cdot(-1)^k\binom{j}{k}\\
&= \sum_{k=0}^{n-1}(-1)^{i+k}\binom{j}{i}\binom{j-i}{k-i}\\
&= (-1)^i\binom{j}{i}\sum_{k=0}^{n=1}(-1)^k\binom{j-i}{k-i}\\
&= (-1)^i\binom{j}{i}(1-1)^{j-i}\\
&= [i=j],
\end{aligned}
\]
所以 \(V=V^{-1}\)。因此答案为
\[V\begin{bmatrix}
\lambda_0^m\\
&\lambda_1^m\\
&&\ddots\\
&&&\lambda_n^m
\end{bmatrix}V\boldsymbol p.
\]
其中 \(V\) 的变换效果是一个差卷积，NTT 实现即可。复杂度 \(\mathcal O(n\log n)\)。
\(\mathcal{Code}\)
/*+Rainybunny+*/
#include <bits/stdc++.h>
#define rep(i, l, r) for (int i = l, rep##i = r; i <= rep##i; ++i)
#define per(i, r, l) for (int i = r, per##i = l; i >= per##i; --i)
typedef long long LL;
const int MAXL = 1 << 18, MOD = 998244353, G = 3;
int n, p[MAXL + 5], fac[MAXL + 5], ifac[MAXL + 5];
LL m;
inline int sgn(const int u) { return u & 1 ? MOD - 1 : 1; }
inline int mul(const int u, const int v) { return 1ll * u * v % MOD; }
inline int sub(int u, const int v) { return (u -= v) < 0 ? u + MOD : u; }
inline int add(int u, const int v) { return (u += v) < MOD ? u : u - MOD; }
inline int mpow(int u, int v) {
int ret = 1;
for (; v; u = mul(u, u), v >>= 1) ret = mul(ret, v & 1 ? u : 1);
return ret;
}
inline void init() {
fac[0] = 1;
rep (i, 1, n) fac[i] = mul(i, fac[i - 1]);
ifac[n] = mpow(fac[n], MOD - 2);
per (i, n - 1, 0) ifac[i] = mul(i + 1, ifac[i + 1]);
}
inline void ntt(const int n, int* u, const int tp) {
static int rev[MAXL + 5]; int lgn = 31 - __builtin_clz(n);
rep (i, 1, n - 1) rev[i] = rev[i >> 1] >> 1 | (i & 1) << lgn >> 1;
rep (i, 0, n - 1) if (i < rev[i]) std::swap(u[i], u[rev[i]]);
for (int stp = 1; stp < n; stp <<= 1) {
int wi = mpow(G, (MOD - 1) / (stp << 1));
for (int j = 0; j < n; j += stp <<1 ) {
for (int wk = 1, k = j; k < j + stp; ++k, wk = mul(wk, wi)) {
int ev = u[k], ov = mul(wk, u[k + stp]);
u[k] = add(ev, ov), u[k + stp] = sub(ev, ov);
}
}
}
if (!~tp) {
std::reverse(u + 1, u + n);
int inv = mpow(n, MOD - 2);
rep (i, 0, n - 1) u[i] = mul(u[i], inv);
}
}
inline void transV() {
static int f[MAXL + 5], g[MAXL + 5];
int len = 1 << 32 - __builtin_clz((n << 1) - 1);
rep (i, 0, len - 1) f[i] = g[i] = 0;
rep (i, 0, n - 1) f[i] = mul(fac[i], p[i]), g[i] = ifac[n - 1 - i];
ntt(len, f, 1), ntt(len, g, 1);
rep (i, 0, len - 1) f[i] = mul(f[i], g[i]);
ntt(len, f, -1);
rep (i, 0, n - 1) p[i] = mul(mul(sgn(i), ifac[i]), f[n - 1 + i]);
}
int main() {
scanf("%d %lld", &n, &m);
++n, m %= MOD - 1, init();
rep (i, 0, n - 1) scanf("%d", &p[i]);
transV();
rep (i, 0, n - 1) p[i] = mul(p[i], mpow(i + 1, MOD - 1 - m));
transV();
rep (i, 0, n - 1) printf("%d%c", p[i], i < repi ? ' ' : '\n');
return 0;
}