官方题解只能做到 $O(n^5)$？玩的太差了。

你的 GPT 只能做到 $O(n^3)$？玩的太差了。

下面给出一个 $O(n^2\log^2 n)$ 时间、$O(n^2)$ 空间的算法；它与输出规模的下界 $\Omega(n^2)$ 只差对数因子，在忽略对数因子的意义下已经最优。([QOJ][1])

关键恒等式与代数化简

设 $H_\ell(x,y)=\sum_{p\in S_\ell}x^{\operatorname{exc}(p)}y^{\operatorname{des}(p)}$，其中 $\operatorname{exc}$ 表示超过数，$\operatorname{des}$ 表示降低数，并规定 $H_0(x,y)=1$。题目要求的就是 $H_n(x,y)$ 的所有系数。

这里使用 Foata–Han 关于这两个统计量联合分布的生成函数恒等式，即其论文 Fix-Mahonian Calculus III; a Quadruple Distribution 中的式 $(1.15)$。注意，它确实描述的是“超过数与降低数”，而不是“排列与逆排列的降低数”。改用我们的记号，恒等式为：([arXiv][2])

$$ \sum_{\ell\ge0}\frac{H_\ell(x,y)}{(1-y)^{\ell+1}}z^\ell = \sum_{r\ge0}y^r \frac{(1-x)(1-xz)^r} {(1-z)\bigl((1-z)^r-x(1-xz)^r\bigr)}. $$

以下从这个恒等式推导算法。

将右侧 $y^r$ 对应的有理函数的 $z^n$ 系数记为 $C_{n,r}(x)$，于是 $H_n(x,y)=(1-y)^{n+1}\sum_{r\ge0}C_{n,r}(x)y^r$。这相当于对降低数这一维做了一次二项式变换。因为最后只需要 $y^0,\ldots,y^{n-1}$ 的系数，只用计算 $C_{n,0}(x),\ldots,C_{n,n-1}(x)$。

直接逐个展开这些有理函数仍然不够快。关键是把 $C_{n,r}(x)$ 写成关于 $r$ 的多项式，并一次性算出它的系数。

令 $R=(1-xz)/(1-z)$。上面的有理函数可以展开成 $\frac{1-x}{1-z}\sum_{a\ge0}x^aR^{r(a+1)}$。又因为 $R=1+\frac{(1-x)z}{1-z}$，使用二项式定理，再利用 $[z^{n-j}](1-z)^{-j-1}=\binom nj$，得到

$$ C_{n,r}(x) = \sum_{j=0}^{n} \binom nj(1-x)^{j+1} \sum_{a\ge0}x^a\binom{r(a+1)}j. $$

现在引入有符号第一类斯特林数 $c(j,i)$，定义为 $t(t-1)\cdots(t-j+1)=\sum_{i=0}^jc(j,i)t^i$。它满足 $c(0,0)=1$，以及 $c(j,i)=c(j-1,i-1)-(j-1)c(j-1,i)$。因此 $\binom{r(a+1)}j=\frac1{j!}\sum_{i=0}^jc(j,i)r^i(a+1)^i$。

另一方面，普通欧拉多项式 $A_i(x)$ 满足 $\sum_{a\ge0}(a+1)^ix^a=A_i(x)/(1-x)^{i+1}$，其中 $A_0(x)=1$；对于 $i\ge1$，它就是按降低数计数的 $i$ 阶排列生成多项式。([arXiv][2])

把这两个展开式代入并交换求和次序，得到 $C_{n,r}(x)=\sum_{i=0}^nr^iA_i(x)\sum_{j=i}^n\binom nj\frac{c(j,i)}{j!}(1-x)^{j-i}$。

定义 $B_i(z)=\sum_{j=i}^n\binom nj\frac{c(j,i)}{j!}z^{j-i}$，再定义 $D_i(x)=A_i(x)B_i(1-x)$，整个式子就变成了

$$ \boxed{C_{n,r}(x)=\sum_{i=0}^{n}r^iD_i(x).} $$

这就是算法最重要的转化：先通过多项式平移和乘法构造所有 $D_i$，随后把系数矩阵转置，通过多点求值得到所有 $C_{n,r}$。 不需要进行三维动态规划。

这里还有一个有用的次数界。对于 $i\ge1$，$\deg A_i=i-1$，而 $\deg B_i\le n-i$，所以 $\deg D_i\le n-1$；对于 $i=0$，由 $c(j,0)=0$（$j>0$）可知 $D_0(x)=1$。因此所有 $D_i$ 总共只有 $O(n^2)$ 个系数。

如何实现准二次复杂度

记 $\mu(d)$ 为两个长度 $O(d)$ 的多项式相乘的时间复杂度，以免与题目中的模数 $M$ 混淆。

首先在 $O(n^2)$ 时间内预处理所有 $c(j,i)$、阶乘、阶乘逆元，以及普通欧拉多项式 $A_0,\ldots,A_n$。欧拉多项式可以用递推 $A_{i+1}(x)=(1+ix)A_i(x)+x(1-x)A_i'(x)$ 计算，每个系数的转移都是常数时间。

接下来构造全部 $B_i$。它们的每个系数都是 $\binom njc(j,i)/j!$，直接查预处理表即可，总时间仍然是 $O(n^2)$。

计算 $D_i(x)=A_i(x)B_i(1-x)$ 时，不能逐项展开 $B_i(1-x)$，而要使用快速多项式平移。具体来说，对于 $B(z)=\sum_{j=0}^db_jz^j$，有 $[x^k]B(1-x)=\frac{(-1)^k}{k!}\sum_{j=k}^db_j\frac{j!}{(j-k)!}$。令 $u_{d-j}=b_jj!$、$v_j=1/j!$，计算一次卷积 $w=u*v$，就得到 $[x^k]B(1-x)=(-1)^kw_{d-k}/k!$。

因此，每个 $B_i(1-x)$ 可以通过一次卷积算出，再与 $A_i(x)$ 做一次乘法，得到 $D_i(x)$。全部 $D_i$ 的计算时间为 $O(n\mu(n))$。

然后令 $d_{i,m}=[x^m]D_i(x)$。对于每个固定的超过数 $m$，构造关于新变量 $t$ 的多项式 $F_m(t)=\sum_{i=0}^nd_{i,m}t^i$。根据刚才的核心等式，有 $F_m(r)=[x^m]C_{n,r}(x)$。于是问题转化为：

对每个 $m$，计算一个次数不超过 $n$ 的多项式在 $0,1,\ldots,n-1$ 这 $n$ 个点上的值。

这一步不能逐点使用 Horner 法，否则总复杂度会退化为 $O(n^3)$。使用乘积树与余式树进行快速多点求值：先建立叶子多项式 $t-r$ 的乘积树，再将 $F_m$ 对根多项式 $\prod_{r=0}^{n-1}(t-r)$ 取余，随后沿树递归取余。到达叶子 $t-r$ 时，余数就是 $F_m(r)$。单个多项式的求值复杂度为 $O(\mu(n)\log n)$，全部 $n$ 个多项式共 $O(n\mu(n)\log n)$；所有求值可以共用同一棵乘积树。([SPECFUN][3])

此时已经得到 $c_{r,m}=[x^m]C_{n,r}(x)$，但它还不是答案。最后还要乘回 $(1-y)^{n+1}$，即

$$ \boxed{ h(n,m,k) = \sum_{r=0}^{k} (-1)^{k-r}\binom{n+1}{k-r}\,c_{r,m}. } $$

对于固定的 $m$，这就是序列 $c_{0,m},\ldots,c_{n-1,m}$ 与公共序列 $b_j=(-1)^j\binom{n+1}{j}$ 的一次卷积，只保留前 $n$ 项。全部答案的恢复时间为 $O(n\mu(n))$。

例如，当 $n=3$ 时，算出的三个辅助多项式为 $C_{3,0}=1$、$C_{3,1}=4+3x+x^2$、$C_{3,2}=10+13x+4x^2$。乘回 $(1-y)^4$ 并保留 $y$ 的前 $3$ 项，就得到 $H_3(x,y)=1+(3x+x^2)y+xy^2$，与样例一致。

综上，算法总时间为 $O(n^2+n\mu(n)\log n)$。在使用 $\mu(n)=O(n\log n)$ 的快速卷积时，就是 $O(n^2\log^2 n)$；系数表占用 $O(n^2)$ 空间，乘积树占用 $O(n\log n)$ 空间，总空间为 $O(n^2)$。题目要求显式输出 $n^2$ 个整数，因此有 $\Omega(n^2)$ 的输出下界，这个算法达到了忽略对数因子后的最优复杂度。

实现时，题目保证 $n\le60$，且 $M$ 是 $10^8$ 到 $10^9$ 之间的质数，因此涉及的阶乘均可求逆，但不能假设 $M$ 本身适合 NTT。([QOJ][4]) 附带实现使用三个辅助 NTT 素数，通过中国剩余定理恢复整数卷积后再对 $M$ 取模；在题目范围内，三个素数的乘积大于卷积系数的上界，所以这里是精确运算，不涉及浮点误差。