人类只能做到 $O(n^6)$？玩的太差了！可以做到 $O(n^4)$ 时间、$O(n^2)$ 额外空间，不需要指数状态、NTT 或多项式插值。关键是固定最终的 $n$，利用所有目标多项式满足同一个递推的性质，直接求出整个解空间，而不是不断递推更小规模的问题。
先说明复杂度结论的边界：题目需要输出 $\Theta(n^3)$ 个数，因此存在 $\Omega(n^3)$ 的输出下界。下面完整证明并实现的是 $O(n^4)$ 算法，不将它声称为已经证明的渐近最优算法。
从排列到多项式
记 $T_k=k(k+1)/2$，$N=T_n$。题目中的对象是所有单元素集合与双元素集合；对于 $i< j$，它们必须满足 ${i}<{i,j}<{j}$，并且所有单元素集合依次出现。数据范围为 $n\le40$，模数是 $10^8< p< 10^9$ 的质数。([QOJ][1])
给每个对象分配一个位于 $(0,1)$ 的实数，按照实数从小到大排序。每个固定排列对应的区域体积都是 $1/N!$。
设单元素集合 ${i}$ 对应的实数为 $t_i$。固定 $0< t_1 < \cdots< t_n< 1$ 后，${i,j}$ 对应的实数可以独立地取在 $(t_i,t_j)$ 中。因此，把所有双元素集合对应的变量积分掉，剩下的权重恰好是范德蒙德乘积 $\Delta(t)=\prod_{i< j}(t_j-t_i)$。
设 $Z=\int_{0< t_1 <\cdots< t_n< 1}\Delta(t),dt$，那么总合法排列数就是 $N!Z$。
对 $0\le k\le n$，定义区域 $D_k(x)={0< t_1< \cdots< t_k< x< t_{k+1}<\cdots< t_n< 1}$，并令 $E_k(x)=Z^{-1}\int_{D_k(x)}\Delta(t),dt$。
再设 $B_{k,d}$ 表示“前 $d$ 项中恰好有 $k$ 个单元素集合”的合法排列数。一个固定排列在恰有前 $d$ 个实数小于 $x$ 时，对应的体积为 $x^d(1-x)^{N-d}/(d!(N-d)!)$。于是 $E_k(x)=Z^{-1}\sum_{d=0}^N B_{k,d}x^d(1-x)^{N-d}/(d!(N-d)!)$。
直接处理这种基底不方便。令 $x=y/(1+y)$，定义 $P_k(y)=(1+y)^N E_k(y/(1+y))$，就得到：
$$ B_{k,d}=Z\,d!\,(N-d)!\,[y^d]P_k(y). $$
我们最终要求的是 ${i}$ 出现在第 $j$ 项的方案数。令 $Q_i(y)=\sum_{k=0}^{i-1}P_k(y)$，再令 $C_{i,d}=Z,d!,(N-d)!,[y^d]Q_i(y)$。那么 $C_{i,d}$ 正是“${i}$ 尚未出现在前 $d$ 项”的合法排列数，因此答案为 $A_{i,j}=C_{i,j-1}-C_{i,j}$。
下面只需要计算这些 $P_k$。
它们有非常重要的首尾次数限制。前缀中恰好有 $k$ 个单元素集合时，所有完全包含于 $[k]$ 的集合都已经出现，因此 $d\ge T_k$。类似地，后缀至少包含 $T_{n-k}$ 个对象。记 $b_k=N-T_{n-k}$，则 $P_k$ 只有 $[T_k,b_k]$ 内的系数可能非零。
将排列反转，同时把编号 $i$ 换成 $n+1-i$，又得到 $P_k(y)=y^NP_{n-k}(1/y)$。
还需要知道首项系数。定义
$$ S_m(a)=\int_{0< u_1 <\cdots < u_m< 1}\Delta(u)\prod_{i=1}^m u_i^a\,du =\frac{1}{\prod_{i=1}^m(a+i)} \prod_{1\le i< j\le m}\frac{j-i}{2a+i+j}. $$
这里 $S_0(a)=1$。这个显式乘积是 Selberg 积分公式在相应参数下的特例，已经完全化成小整数乘除，不需要在程序中实现任何积分操作。([arXiv][2])
显然 $Z=S_n(0)$。当 $x\to0$ 时，将前 $k$ 个变量写成 $t_i=xu_i$，它们自身的积分贡献 $x^{T_k}S_k(0)$；跨越两组的差值则让后面的每个变量多出一个 $k$ 次方因子。因此，记 $\alpha_k=[y^{T_k}]P_k(y)$，有 $\alpha_k=S_k(0)S_{n-k}(k)/Z$。由反转对称性，$[y^{b_k}]P_k(y)=\alpha_{n-k}$。
这些乘积里的分子、分母都是非零的小整数，在题目模数下全部可逆，特别地，每个 $\alpha_k$ 都非零。
固定 $n$，直接求出整个解空间
这里引入的辅助量数量只有 $n+1$。
固定一个 $k$，暂时省略下标 $k$。令 $e_r$ 为 $n$ 个变量的 $r$ 次基本对称多项式，即从这些变量中选出 $r$ 个相乘后求和，且 $e_0=1$。定义 $J_r(x)=\bigl(Z\binom nr\bigr)^{-1}\int_{D_k(x)}\Delta(t)e_r(x-t_1,\ldots,x-t_n),dt$。
特别地，$J_0=E_k$。
记 $a_r=(n-r)(n-r+1)$、$c_r=(n-r)(2n-r+2)$、$d_r=r(n-r+2)$，以及 $\sigma=x(1-x)$。这些辅助积分满足：
$$ c_rJ_{r+1} =a_r(2x-1)J_r+2\sigma J_r'-d_r\sigma J_{r-1}, \qquad 0\le r\le n, $$
其中约定 $J_{-1}=J_{n+1}=0$。这是 Selberg 积分微分递推在本题参数下的形式。([arXiv][2])
这个关系也可以直接用分部积分验证。设 $e_r^{(i)}$ 表示删去变量 $x-t_i$ 后的基本对称多项式，对 $\sum_i\partial_{t_i}\bigl(t_i(1-t_i)\Delta(t)e_r^{(i)}\bigr)$ 在 $D_k(x)$ 上积分。端点 $0,1$ 的边界项被 $t_i(1-t_i)$ 消去，两个变量重合处的边界项被 $\Delta$ 消去；剩下的移动边界正好对应 $J_r'$ 中的边界贡献。展开时利用 $\sum_i e_r^{(i)}=(n-r)e_r$、$\sum_i(x-t_i)e_r^{(i)}=(r+1)e_{r+1}$，并将 $\Delta$ 的导数按变量对配对，即得到上式。对于 $r=n$，关系退化为 $J_n'=nJ_{n-1}$。
关键在于：递推的系数与 $k$ 完全无关。 不同的 $P_k$ 只是同一个线性系统的不同解。
令 $H_r(y)=(1+y)^{N+r}J_r(y/(1+y))$，并记 $h_{r,\ell}=[y^\ell]H_r(y)$。注意 $H_0=P_k$。代入上面的微分递推并比较系数，得到：
$$ c_rh_{r+1,\ell} =(2\ell-a_r)h_{r,\ell} +\bigl(2\ell-2-r(2n-r+3)\bigr)h_{r,\ell-1} -d_rh_{r-1,\ell-1}. $$
所有次数为负的系数，以及下标越界的辅助量，都视为零。
这个式子只涉及当前次数和上一个次数，因此可以按照 $\ell=0,1,\ldots,N$ 依次计算。
如果 $\ell$ 不是三角数，那么所有 $2\ell-a_r$ 都非零。此时从 $r=n$ 向 $r=0$ 倒推，就能算出这一次数的全部辅助量。
如果 $\ell=T_s$，则恰好在 $r_*=n-s$ 处有 $2\ell-a_{r_*}=0$。这时，将 $h_{0,\ell}$ 作为一个自由参数：从 $r=0$ 开始正向递推，算到 $h_{r_*,\ell}$；另外从 $r=n$ 开始反向递推，算到 $h_{r_*+1,\ell}$。这样就避开了唯一的零分母。
为什么这里可以自由指定 $h_{0,T_s}$，而被跳过的方程不会产生额外限制？
因为我们已经知道，真正的 $P_0,\ldots,P_n$ 都满足这个系统，并且它们的最低次数分别是 $T_0,\ldots,T_n$，对应系数分别为非零的 $\alpha_0,\ldots,\alpha_n$。因此，它们在这 $n+1$ 个次数上的取值矩阵是对角线非零的三角矩阵。
也就是说，任意指定 $[y^{T_0}]H_0,\ldots,[y^{T_n}]H_0$，都恰好对应某个 $P_0,\ldots,P_n$ 的线性组合。上述递推在普通次数处唯一，在三角数次数处恰有一个自由度，所以生成的正是这个线性组合。被跳过的方程必然成立，而不是通过随意赋值掩盖零分母。
于是我们得到了一个生成器：输入一个多项式在全部三角数次数上的 $n+1$ 个系数，即可在 $O(nN)=O(n^3)$ 时间内依次生成它的全部系数。
接下来恢复真正的 $P_k$。
先用生成器构造一组基底 $g_0,\ldots,g_n$，满足 $[y^{T_t}]g_s=[s=t]$。不需要保存所有基底的全部系数，只保存矩阵 $R_{s,l}=[y^{b_l}]g_s$，它只有 $O(n^2)$ 个元素。
由于 $g_k$ 在 $T_k$ 之前的系数都是零，它一定属于 $P_k,\ldots,P_n$ 的线性张成空间。将其乘上 $\alpha_k$，就使最低次项与 $P_k$ 相同。
现在按照 $k=n,n-1,\ldots,0$ 处理。对于当前多项式，依次用已经恢复出的 $P_n,P_{n-1},\ldots,P_{k+1}$ 消去次数 $b_n,b_{n-1},\ldots,b_{k+1}$ 处的系数。消去 $b_l$ 处系数时，除数就是已知非零的 $\alpha_{n-l}$。
这些 $P_l$ 的最高次数严格递增，因此这相当于三角消元。所有 $l>k$ 对应的成分都消去后，剩下的只能是 $P_k$；而最低次项已经正确，所以比例系数也是正确的。
消元时同步维护变换矩阵 $U$，使 $P_k=\sum_sU_{k,s}g_s$。由于 $g_s$ 的定义，实际上 $U_{k,s}=[y^{T_s}]P_k$，也就是重新生成 $P_k$ 所需的自由参数。
最后，为了按题目要求逐行输出，又不保存整张答案表，对每个 $i$ 把 $U$ 的前 $i$ 行相加，得到 $Q_i=\sum_{k< i}P_k$ 的自由参数，再调用一次生成器。每生成一个次数 $d$ 的系数，就计算 $C_{i,d}$，与上一个值作差并直接输出。
构造全部基底需要 $O(n^2N)$ 时间；小矩阵消元需要 $O(n^3)$ 时间；重新生成全部输出行需要 $O(n^2N)$ 时间。总时间为 $O(n^4)$。生成器只保留相邻两个次数的辅助量，其他矩阵和阶乘数组也都只有 $O(n^2)$ 大小，因此额外空间为 $O(n^2)$。
所有需要求逆的小整数的绝对值都不超过 $4N$，可以线性预处理逆元，不需要对每次除法做快速幂。