官方题解只能做到 $O(n \cdot k^{1.5})$？玩的太差了。这道题可以先做到 $O(n^2k)$，再利用“只需要最后一项”这一点，对固定的 $n$ 进一步做到 $O(\sqrt{k}\log k)$。关键是把最大值的统计转化成只有 $n(n+1)/2$ 个状态的递推，再对这个递推进行批量跳转。
题目要求计算 $n^k \mathbb E[ \max(a_1, \ldots, a_n)]$，其中 $n\le 10$、$k\le 10^5$，模数为 $p=998244353$。由于每个操作序列的概率都是 $n^{-k}$，要求的实际上就是所有操作序列的最大值之和。
下面先在有理数形式幂级数中推导，最后将需要的低阶系数模 $p$。计算过程中涉及的整数分母均小于 $p$。$n=1$ 时答案直接为 $k$，以下考虑 $n\ge 2$。
把最大值统计压缩成一个小递推
记进行 $t$ 次操作时的答案为 $A_t$，并定义指数生成函数 $A(x)=\sum_{t\ge 0}A_tx^t/t!$。再记 $u_j(x)=x^j/j!$、$E_j(x)=\sum_{i=0}^ju_i(x)$，约定 $E_{-1}(x)=0$。
对于一个最终状态 $(a_1,\ldots,a_n)$，产生它的操作序列数为 $t!/(a_1!\cdots a_n!)$。因此，$E_j(x)^n-E_{j-1}(x)^n$ 正好统计最大值等于 $j$ 的状态，得到 $A(x)=\sum_{j\ge 0}j\bigl(E_j(x)^n-E_{j-1}(x)^n\bigr)$。
直接对每个 $j$ 计算这个表达式并不好。我们改为记录“有多少个变量同时取得最大值”，并引入最大值的少量幂次。
定义 $H_{r,s}(x)=\sum_{j\ge 0}j^s u_j(x)^rE_{j-1}(x)^{n-r}$，其中 $1\le r\le n$，约定 $0^0=1$。这里选定的 $r$ 个变量都等于 $j$，其他变量严格小于 $j$，额外权重为 $j^s$。乘上 $\binom nr$ 后，才是任选这 $r$ 个变量的统计。
展开 $E_j^n-E_{j-1}^n$，立即得到 $A(x)=\sum_{r=1}^n\binom nrH_{r,1}(x)$。
乍看之下，$s$ 可能需要无限增长。真正有用的是下面这个闭合关系：只需要保留 $0\le s< r$ 的状态。
令 $G_{r,s}(x)=x^{-s}H_{r,s}(x)$。当 $0\le s\le r$ 时，这仍然是普通形式幂级数：对于 $s>0$，$j=0$ 项消失，其余项至少含有 $x^r$。
注意到 $j^ru_j^r=x^ru_{j-1}^r$。将 $H_{r,r}$ 中的下标减一，再使用 $E_j=E_{j-1}+u_j$ 展开，就得到
$$ G_{r,r}(x)=\sum_{a=r}^n\binom{n-r}{a-r}G_{a,0}(x). $$
也就是说，每当第二维增长到第一维时，都能用第二维为 $0$ 的状态消去。因此，需要真正存储的状态只有 $(r,s)$，其中 $1\le r\le n$、$0\le s < r$，总数为 $D=n(n+1)/2$。
接下来求这些状态之间的递推。
记 $\theta=x\frac{\mathrm d}{\mathrm dx}$。由 $\theta u_j=ju_j$ 以及 $\theta E_{j-1}=xE_{j-1}-ju_j$，对 $H_{r,s}$ 求导可得 $\theta H_{r,s}=rH_{r,s+1}+(n-r)xH_{r,s}-(n-r)H_{r+1,s+1}$。换成 $G$ 后就是 $(\theta+s)G_{r,s}=x\bigl(rG_{r,s+1}+(n-r)G_{r,s}-(n-r)G_{r+1,s+1}\bigr)$。
这已经是一个只含 $D$ 个未知函数的封闭线性微分系统。
为了把最终答案所需的阶乘也吸收进递推，定义 $w_{r,s}(t)=(t+1)![x^t]G_{r,s}(x)$。比较系数，得到
$$ w_{r,s}(t)= \frac{t+1}{t+s} \left( r\,w_{r,s+1}(t-1) +(n-r)w_{r,s}(t-1) -(n-r)w_{r+1,s+1}(t-1) \right). $$
这里 $t\ge 1$。当 $s=r-1$ 时，右侧的 $w_{r,r}$ 使用前面的闭合关系计算；当 $r=n$ 时，最后一项直接省略，不访问不存在的状态。
初值也非常简单：$w_{n,0}(0)=1$，其余存储状态均为 $0$。 原因是常数项只能来自所有变量都为 $0$ 的情况。
每一轮先计算所有辅助值 $w_{r,r}=\sum_{a=r}^n\binom{n-r}{a-r}w_{a,0}$，再用上一轮状态计算下一轮。辅助值总共需要 $O(n^2)$ 次运算，真正的状态转移也只需要 $O(n^2)$ 次运算。所有状态都只依赖上一轮，没有同轮依赖问题。
最后考虑如何读出答案。对于 $r\ge 2$，有 $H_{r,1}=xG_{r,1}$；对于 $r=1$，则使用闭合关系 $H_{1,1}=xG_{1,1}=x\sum_{a=1}^n\binom{n-1}{a-1}G_{a,0}$。再利用 $n\binom{n-1}{r-1}=r\binom nr$，得到
$$ A_k= \sum_{r=1}^n\binom nr \left( r\,w_{r,0}(k-1) +[r\ge 2]\,w_{r,1}(k-1) \right). $$
因此，直接递推的时间复杂度为 $O(n^2k)$。滚动数组只占 $O(n^2)$ 空间；使用线性预处理的逆元表时，总空间为 $O(k+n^2)$。这里最大的分母是 $k+n-2$，不会出现模意义下不可逆的情况。
但这并不是关于 $k$ 的最终复杂度：我们没有必要把中间的 $k-1$ 轮全部算出来。
只计算最后一项：等距插值与矩阵分块
把所有 $w_{r,s}(t)$ 排成列向量 $W_t$。前面递推中括号里的线性组合可以写成 $BW_{t-1}$，其中 $B$ 是一个与 $t$ 无关的常数矩阵。辅助状态 $w_{r,r}$ 已通过闭合关系展开，所以 $B$ 的维数就是 $D$，且只有 $O(n^2)$ 个非零元素。
再令 $J$ 为对角矩阵，状态 $(r,s)$ 对应的对角元为 $s$，则递推可以写成 $W_t=(t+1)(tI+J)^{-1}BW_{t-1}$。
这是一个矩阵乘积问题，不过转移矩阵的元素是有理函数。先把分母统一。
定义 $q(z)=\prod_{\substack{0\le s< n\\s\ne 1}}(z+s)$，并令 $P(z)=q(z)(z+1)(zI+J)^{-1}B$。那么 $P(z)$ 是次数至多 $d=n-1$ 的多项式矩阵。对于 $s=1$ 的行，分母 $z+1$ 被分子抵消；对于其他行，分母 $z+s$ 被 $q(z)$ 抵消。
于是 $W_t=P(t)W_{t-1}/q(t)$。
为了同时计算分子和分母，可以再增加一维，使用分块对角矩阵 $\widehat P(z)=\operatorname{diag}(P(z),q(z))$。从向量 $(W_0,1)$ 出发，计算 $\widehat P(K)\cdots\widehat P(1)$ 的作用，其中 $K=k-1$。最后一维给出总分母，前 $D$ 维除以它，就得到 $W_K$。
至此，问题变成：求一个固定大小、固定次数的多项式矩阵在连续整数点上的有序乘积。
这一类乘积可以用 baby-step/giant-step 加等距采样点平移处理；相比先求出整个矩阵多项式再做普通多点求值，直接维护采样值可以省去额外的对数因子。Bostan–Gaudry–Schost 的多项式系数递推算法使用的正是这一思路。([SPECFUN][2])
下面把本题需要的算法具体写出来。
选择一个为 $2$ 的幂的块长 $b$，定义 $F_m(z)=\widehat P(z+m-1)\cdots\widehat P(z)$。它的次数至多为 $dm$，并且满足 $F_{2m}(z)=F_m(z+m)F_m(z)$。这里的矩阵乘法顺序不能交换。
从 $m=1$ 开始倍增，但不存储 $F_m$ 的系数，而是存储它在 $1,1+b,\ldots,1+dmb$ 上的值。令 $L=dm$，也就是维护 $F_m(1+bj)$，其中 $0\le j\le L$。
为了倍增到 $2m$，只需获得两组值：$F_m(1+bj)$ 和 $F_m(1+m+bj)$，其中 $0\le j\le 2L$。将 $F_m(1+bz)$ 看成关于 $z$ 的多项式，这分别是把已有采样延长，以及把采样点平移 $m/b$。
两组值准备好后，逐点做一次矩阵乘法，就得到了 $F_{2m}(1+bj)$。
这里用到的采样平移可以由拉格朗日插值直接实现。对于次数至多为 $L$ 的多项式，已知 $f(0),\ldots,f(L)$，有 $f(y)=\prod_{j=0}^L(y-j)\sum_{i=0}^L\frac{(-1)^{L-i}f(i)}{i!(L-i)!(y-i)}$。
取 $y=\alpha+t$ 后，求和部分是序列 $\frac{(-1)^{L-i}f(i)}{i!(L-i)!}$ 与序列 $\frac1{\alpha+t-i}$ 的卷积；外面的连乘也可以线性递推。因此，一次卷积便能求出 $O(L)$ 个平移后的采样值，时间为 $O(\mathsf M(L))$，其中 $\mathsf M(L)$ 表示多项式乘法复杂度。遇到目标点恰好是已有采样点时直接复制，不能对零分母求逆。
对矩阵的每个元素应用这一操作，再逐点相乘。如果所采用的矩阵乘法复杂度为 $O(D^\omega)$，一次倍增的代价就是 $O(D^\omega L+D^2\mathsf M(L))$。各层的 $L$ 按倍数增长，所以构造长度为 $b$ 的块，总共只需要 $O(D^\omega db+D^2\mathsf M(db))$。
得到 $F_b$ 后，第 $j$ 个完整块的转移就是 $F_b(1+jb)$。如果完整块的数量超过已有采样数，就继续用同一个采样平移算法，每次补出 $O(db)$ 个块转移。
随后，按照块的先后顺序，将这些矩阵作用到当前向量上即可。这里只需要矩阵乘向量，不需要把所有块矩阵先乘成一个总矩阵。 不足一块的末尾部分直接处理。
取 $\mathsf M(L)=O(L\log L)$，整个算法的时间可以写成
$$ T(b)=O\left( D^\omega db+ D^2\left(db+\frac Kb\right)\log(db+2) \right). $$
这也说明，块长不一定应该机械地取 $\sqrt K$：构造块时做的是矩阵乘法，应用块时做的是矩阵乘向量，两者的代价不同。
忽略对数因子，选择 $b=\Theta\left(\sqrt{K/(dD^{\omega-2})}\right)$，并取附近的 $2$ 的幂，可以得到 $\widetilde O\left(D^{(\omega+2)/2}\sqrt{dK}\right)$。若这个块长小于 $1$，直接使用前面的稀疏递推即可。
代入 $D=O(n^2)$、$d=O(n)$，得到 $\widetilde O(n^{\omega+5/2}\sqrt k)$ 的确定性算法。使用普通矩阵乘法，即 $\omega=3$ 时，一个明确的上界是 $O(n^{11/2}\sqrt k\log(nk))$；对于题目中的固定小常数 $n$，就是 $O(\sqrt k\log k)$。
额外采样可以分批计算、立即作用到向量后丢弃，所以空间为 $O(D^2db)$，固定 $n$ 时为 $O(\sqrt k)$。这个版本不能再预处理一直到 $k$ 的阶乘或逆元表，否则会重新引入 $O(k)$ 时间；前面把阶乘吸收进 $w_{r,s}(t)$，正好避免了这一点。
综合起来，可以在两种算法中选择代价较小的一个，得到 $O\left(\min\{n^2k,\ n^{11/2}\sqrt k\log(nk)\}\right)$ 的上界。对于 $n=10,k=10^5$ 这样的范围，稀疏递推的常数很小；但在固定 $n$、考察 $k$ 的渐近增长时，分块插值确实把线性时间降到了次线性时间。
提交记录：https://qoj.ac/submission/2936404 轻松比人类快 100 倍。