官方题解只能做到 $O(n^3)$？玩的太差了。这题可以做到 $O(n^2\log n)$ 时间、$O(n^2)$ 空间，不必停在立方复杂度的计数 DP。关键是先把“不同树可能产生相同 DFS 序”的去重问题转化成一个形式幂级数方程，再通过倍增求解；其中每轮求解的矩阵也要用卷积批量构造，否则仍然会退化成 $O(n^3)$。
以下把输入模数记作 $p$。题目保证 $n\le 800$，且 $p$ 是大于 $n$ 的素数，因此后面出现的整数、阶乘的除法都可以用逆元实现。([QOJ][1])
从不同 DFS 序到计数方程
先对原树进行“左儿子、右兄弟”转换：每个点的最小儿子作为左儿子，下一个兄弟作为右儿子。原来的父子关系满足父亲编号更小，兄弟又按编号递增排列，所以转换后，二叉树的每条边都从较小编号指向较大编号。
原来的 DFS 序恰好就是这棵二叉树的先序遍历。删去只起辅助作用的根 $1$，问题便等价于：
统计有多少个长度为 $n-1$ 的排列，能够成为某棵递增二叉树的先序遍历。
反方向的转换也成立，因此这里没有遗漏或增加排列。这也是文献中 sorted permutations 的一个等价计数模型。([ScienceDirect][2])
接下来直接对排列计数，不再对树计数。
对于一个由不同数字组成的序列，可以贪心尝试将它还原成递增二叉树。维护一个递增栈，初始放入第一个数字。处理下一个数字 $x$ 时：
- 若 $x$ 大于栈顶，就将它接成刚访问节点的左儿子，然后入栈。
- 否则，弹出所有大于 $x$ 的节点。如果栈仍非空，再弹出一个节点，将 $x$ 接成它的右儿子，然后把 $x$ 入栈。
- 如果弹出所有大于 $x$ 的节点后栈为空，则当前这棵树已经无法容纳 $x$，只能从 $x$ 开始一棵新树。
这里始终选择最深的可用位置。选择更高的祖先，只会额外关闭若干尚可使用的分支，而不会保留更多后续选择；因此，贪心失败时，其他连接方式也不可能把此前的整个序列放进同一棵树。于是，一个非空排列合法，当且仅当整个过程中不需要另起一棵树。
令 $c(\pi)$ 表示处理排列 $\pi$ 时“另起一棵树”的次数，空排列也定义为 $c(\varnothing)=0$。这个统计量的好处在于，它在最小值处分解时有很简单的递推。
将非空排列写成 $\pi=A,\min,B$。如果 $A$ 非空，读到全局最小值时一定需要另起一棵树，贡献 $c(A)+1$。随后处理 $B$ 时，这个最小值可以提供左、右两个分支：$B$ 的第一棵树放在左边，第一次需要另起树时放在右边，之后的失败才真正增加计数。因此，无论 $A$ 是否为空，都有 $c(\pi)=c(A)+[A\ne\varnothing]+\max(c(B)-1,0)$。
定义指数生成函数 $F(x,y)=\sum_{m\ge0}\frac{x^m}{m!}\sum_{\pi\in S_m}y^{c(\pi)}$，并记 $B(x)=F(x,0)$。于是 $B$ 统计的正是我们需要的合法排列，空排列也计入其中。
左半部分 $A$ 的贡献是 $1+y(F-1)$，右半部分 $B$ 对应统计量减一、但不小于零，其贡献是 $B+(F-B)/y$。在最小值处分解得到
$$ y\frac{\partial F}{\partial x} =(1-y+yF)\bigl(F-(1-y)B\bigr). $$
这里使用指数生成函数，是因为左右两边需要从其余数字中选择各自使用的标号，乘法自动包含相应的二项式系数。
直接展开这个方程并不够快。下面进行一次关键的线性化。
令 $V=(1-y+yF)^{-1}$，代入可得 $yV_x=(1-y)(1+yB)V-1$。再令 $D(x)=\int_0^x B(t),dt$，以及 $H(x,y)=\exp\bigl(x+(y-1)D(x)\bigr)$，计算后恰好有 $(1-y\partial_x)(HV)=H$。
将算子形式地展开为 $(1-y\partial_x)^{-1}=\sum_{a\ge0}y^a\partial_x^a$，然后令 $x=0$。由于 $H(0,y)=V(0,y)=1$，得到 $\sum_{a\ge0}a!y^a[x^a]H(x,y)=1$。
为方便描述，定义线性算子 $\mathcal L(x^ay^b)=a!y^{a+b}$。最终，只需要求解下面这个单变量未知级数 $D(x)$ 的方程：
$$ \mathcal L\!\left(\exp\bigl(x+(y-1)D(x)\bigr)\right)=1, \qquad \mathrm{ans}=n!\,[x^n]D(x). $$
答案里的下标值得注意：$B(x)$ 的第 $m$ 项统计长度为 $m$ 的二叉树先序排列，而原题对应长度 $n-1$；积分以后，恰好变成取 $D$ 的第 $n$ 项。
倍增求解，并把矩阵构造也降到近二次
假设已经求出 $D\bmod x^m$。令 $N=\min(2m,n+1)$，把未知高次项暂时补零，接下来求一个修正量 $\Delta(x)=\sum_{j=m}^{N-1}e_jx^j$。
因为 $\Delta^2$ 至少从 $x^{2m}$ 开始，所以在当前需要的精度内，$\exp((y-1)\Delta)$ 只需要保留常数项和一次项。也就是说，新旧方程之差对 $\Delta$ 是线性的，二次及以上项不会影响 $y^N$ 之前的系数。
记当前的 $H=\exp(x+(y-1)D)$、$K=(y-1)H$，残差为 $R_s=[y^s](\mathcal L(H)-1)$。我们需要解的方程是
$$ R_s+\sum_{j=m}^{s}J_{s,j}e_j=0, \qquad J_{s,j} =\sum_{r=0}^{s-j}(s-r)!\,[x^{s-j-r}y^r]K, \quad m\le s< N. $$
这个矩阵是下三角矩阵，而且对角线特别简单：$J_{s,s}=s![x^0y^0]K=-s!$。
由于 $p>n$，所有对角元均非零，因此可以从低次到高次依次求解，具体就是 $e_s=\bigl(R_s+\sum_{j=m}^{s-1}J_{s,j}e_j\bigr)/s!$。这一部分只需要 $O((N-m)^2)$ 时间，不需要一般的矩阵求逆。
但到这里还没有完成优化：若逐项按定义计算 $J_{s,j}$，每个矩阵元素又枚举一次 $r$，总复杂度依然是立方。 我们必须同时解决两个批量计算问题。
首先，快速构造 $H$ 的系数表。
写成 $H(x,y)=\sum_{r\ge0}H_r(x)y^r$，就有 $H_r(x)=e^{x-D(x)}D(x)^r/r!$。
先计算 $H_0=e^{x-D}$。这一步甚至不需要快速多项式指数：若 $g=x-D$、$h=e^g$，直接用 $h_0=1$ 和 $h_t=\frac1t\sum_{i=1}^t i g_i h_{t-i}$，花 $O(N^2)$ 时间即可。
之后依次计算 $H_r=H_{r-1}D/r$，每次都是一次多项式乘法。
注意，$\mathcal L$ 会把总次数为 $a+b$ 的单项式 $x^ay^b$ 变成 $y^{a+b}$，所以只需要保留 $a+b
设长度为 $N$ 的多项式乘法复杂度为 $M(N)$，那么整个系数表只需要 $O(NM(N)+N^2)$ 时间、$O(N^2)$ 空间。
得到系数表后，所有残差可以在 $O(N^2)$ 时间内求出，因为 $R_s=\sum_{r=0}^{\lfloor s/2\rfloor}(s-r)![x^{s-r}y^r]H$，这里我们只使用 $s\ge m\ge1$ 的残差。
其次，按矩阵的次对角线批量计算 $J$。
固定 $d=s-j$，定义 $Q_d(z)=\sum_{r=0}^{d}[x^{d-r}y^r]K\cdot z^r$，再定义阶乘多项式 $W(z)=\sum_{t=0}^{N-1}t!z^t$。
于是前面的求和恰好就是普通卷积：
$$ J_{s,s-d}=[z^s]\bigl(W(z)Q_d(z)\bigr). $$
也就是说，一次多项式乘法就可以求出一整条次对角线。
这里 $d$ 只需要枚举 $0,\ldots,N-m-1$，共 $O(N)$ 次卷积。构造所有 $Q_d$ 本身只需要 $O(N^2)$ 时间：若记 $h_{a,b}=[x^ay^b]H$，则 $[x^ay^b]K=h_{a,b-1}-h_{a,b}$，越界项视为零，查表即可。
因此，矩阵构造也是 $O(NM(N)+N^2)$，不会留下隐藏的立方瓶颈。
求出 $\Delta$ 后令 $D\leftarrow D+\Delta$，就得到了 $D\bmod x^N$。从 $D=0$、$m=1$ 开始倍增，直到知道第 $n$ 项。
正确性也随之归纳成立：修正量不影响已知的低次项，当前精度内所有非线性误差都消失，而线性系统的对角元均非零，因此每轮恰好确定下一段唯一正确的系数。
复杂度与实现
一次精度为 $N$ 的倍增，构造 $H$、构造矩阵、解三角系统分别需要 $O(NM(N)+N^2)$、$O(NM(N)+N^2)$ 和 $O(N^2)$ 时间。各轮精度按几何级数增长，所以总复杂度是 $O(nM(n)+n^2)$ 时间、$O(n^2)$ 空间。使用 $M(n)=O(n\log n)$ 的快速卷积，即得到 $O(n^2\log n)$ 时间、$O(n^2)$ 空间。
实现时不能直接假设输入模数支持 NTT。题目只保证模数为给定范围内的素数，并没有保证它具有所需阶数的单位根。([QOJ][1])
可以在 $998244353$、$1004535809$、$469762049$ 三个模数下分别卷积，再用 CRT 还原后对 $p$ 取模。每次卷积的输入系数都先表示成 $[0,p-1]$ 中的整数，则单个未取模系数不超过 $801(p-1)^2 < 8.2\times10^{20}$，远小于三个 NTT 素数的乘积，因而可以精确还原，不需要浮点 FFT。