你的 GPT 只能做到 $O(n^{1.82})$？玩的太差了。可以继续降低复杂度。上一版不必使用小步大步法来应用整块算子；改用低阶微分算子的局部性，可以得到更好的界。

设 $k\times k$ 矩阵乘法的复杂度为 $O(k^\omega)$。新的确定性算法复杂度为

$$ \boxed{\widetilde O\!\left(n^{(\omega+1)/2}\right)\text{ 时间，}\quad O(n\log n)\text{ 空间。}} $$

关键改进是：把一个长多项式拆成许多局部余式，让它们共享同一个小矩阵；处理完后，再用 Hermite 插值恢复完整结果。

概率模型与重链分块

原题中，$W_u$ 独立地取 $1,2,3$，概率分别为 $p_{u,1},p_{u,2},p_{u,3}$；我们要求所有给定的 $T_u< T_v$ 同时成立。限制共有 $n-1$ 条，忽略方向后连通，所以它们组成一棵树。([UOJ][1])

先固定所有 $W_u$。忽略重复抽到的卡片后，下一张新卡是 $u$ 的概率，等于 $W_u$ 除以尚未出现的卡片的权值总和。

因此，可以用相互独立的指数时钟 $E_u\sim\operatorname{Exp}(W_u)$ 代替抽卡过程：最早响的时钟按速率比例产生，去掉它之后，指数分布的无记忆性保证后续过程相同。

令 $X_u=e^{-E_u}$，于是 $T_u< T_v$ 等价于 $X_u>X_v$。当 $W_u=j$ 时，$X_u$ 在 $(0,1)$ 上的密度为 $jx^{j-1}$。对权值取平均，得到密度多项式 $q_u(x)=p_{u,1}+2p_{u,2}x+3p_{u,3}x^2$。

任选一个根。定义 $G_u(x)$ 表示父亲的变量取值为 $x$ 时，$u$ 的子树及父子边均满足条件的概率。设 $F_u(x)=q_u(x)\prod_{v\text{ 是 }u\text{ 的儿子}}G_v(x)$。

若边为父亲 $\to u$，则 $G_u(x)=\int_0^xF_u(t),dt$；若边为 $u\to$ 父亲，则 $G_u(x)=\int_x^1F_u(t),dt$。对根也采用第一种定义，答案就是 $G_{\mathrm{root}}(1)$。

所有消息都是多项式，且 $\deg G_u\le3\operatorname{sz}_u$。直接维护系数即可得到平方算法；QOJ 上现有题解给出的复杂度也是 $O(n^2)$。但只加上 NTT 并不能解决长链：每走一步都扫描一次不断增长的多项式，仍然可能花费 $\Theta(n^2)$。([QOJ][2])

按子树大小进行重链剖分。记 $h(u)$ 为重儿子，将轻儿子的贡献提前乘起来，得到 $P_u(x)=q_u(x)\prod_{v\ne h(u)}G_v(x)$。沿重链自底向上，剩下的操作只有两种：$f\mapsto\int_0^xP_u(t)f(t),dt$，或者 $f\mapsto\int_x^1P_u(t)f(t),dt$。

定义操作的重量为 $d_u=\deg P_u+1$。由于 $d_u\le3(1+\sum_{v\text{ 是轻儿子}}\operatorname{sz}_v)$，而每个点只会被 $O(\log n)$ 个轻子树包含，因此 $\sum_u d_u=O(n\log n)$。

考虑一段总重量为 $D$ 的操作。它们的复合可以写成

$$ T[f](x)=\int_0^xU(x,y)f(y)\,dy+\int_0^1V(x,y)f(y)\,dy, \qquad \deg_{\mathrm{total}}U,\deg_{\mathrm{total}}V< D. $$

单次下限积分对应 $U(x,y)=P(y),V=0$；单次上限积分对应 $U(x,y)=-P(y),V(x,y)=P(y)$。继续复合时，交换积分顺序即可验证这个表示保持成立，总次数也满足上述界。

其中 $V$ 保存了反向边产生的边界贡献。后面的局部化只用于 $U$ 部分，不能把 $V$ 也当成局部算子。

核心优化：将积分算子局部化

先只考虑 $U$ 部分。把提高的次数记为 $d$，定义 $S_d(k)=\sum_{b=0}^{d-1}\frac{U_{d-1-b,b}}{k+b+1}$。那么输入单项式 $x^k$ 对输出 $x^{k+d}$ 的贡献就是 $S_d(k)$。

现在做阶乘换基：对 $f(x)=\sum_kf_kx^k$，记 $\widehat f(x)=\sum_k k!f_kx^k$。换基后，相应的转移系数变成 $R_d(k)=\frac{(k+d)!}{k!}S_d(k)$。

关键是，$R_d(z)$ 是次数小于 $d$ 的普通多项式。因为 $R_d(z)=(z+1)\cdots(z+d)\sum_{b=0}^{d-1}\frac{U_{d-1-b,b}}{z+b+1}$，每个分母都被完整消去。

令 $\theta=x\frac{d}{dx}$，有 $\theta x^k=kx^k$。于是，阶乘换基后的变上限积分部分恰好是微分算子 $L=\sum_{d=1}^D x^dR_d(\theta)$。

这个表示有两个重要性质：$L$ 的微分阶数小于 $D$，并且它使多项式次数最多增加 $D$。

上一版在这里继续做小步大步。实际上，微分阶数小于 $D$ 意味着：要知道输出在某点附近的前若干阶展开，只需要知道输入在该点附近多出不到 $D$ 阶的展开。这正好适合批量处理。

设输入 $f$ 有 $m$ 个系数。取一个二次幂 $K$，使 $D\le K< 2D$，再取 $s=O(1+m/D)$ 个不同的点 $\alpha_i$，满足 $sK\ge m+D$。

对每个点计算局部余式 $p_i(x)=\widehat f(x)\bmod(x-\alpha_i)^{2K}$，所以 $\deg p_i< 2K$。由于 $\widehat f-p_i$ 被 $(x-\alpha_i)^{2K}$ 整除，而 $L$ 的微分阶数小于 $D$，有 $L\widehat f-Lp_i$ 被 $(x-\alpha_i)^{2K-D+1}$ 整除，当然也被 $(x-\alpha_i)^K$ 整除。

因此，$r_i(x)=Lp_i(x)\bmod(x-\alpha_i)^K$ 就是我们需要的输出局部余式。又因为 $\deg L\widehat f< m+D\le sK$，这些余式唯一确定整个 $L\widehat f$。

现在所有 $p_i$ 都是次数小于 $2K$ 的多项式，可以共享同一个矩阵。

具体地，建立 $L$ 从次数小于 $2K$ 的多项式到次数小于 $3K$ 的多项式的矩阵 $A$，它的大小为 $3K\times2K$，非零元素为 $A_{k+d,k}=R_d(k)$。

把所有 $p_i$ 的系数排列成一个 $2K\times s$ 的矩阵 $E$，则 $AE$ 的各列就是所有 $Lp_i$。

这里必须注意：$p_i$ 使用的是全局的 $1,x,x^2,\ldots$ 基，而不是各自中心处的 Taylor 基。 因此它们确实可以乘同一个 $A$。直接把不同中心的 Taylor 系数送入同一个矩阵是不正确的。

每 $K$ 列作为一组，把这个长方形乘法拆成常数次 $K\times K$ 矩阵乘法，代价为 $O(K^\omega\lceil s/K\rceil)=O(D^\omega+mD^{\omega-2})$。

余式的求值和恢复都可以使用快速 Hermite 求值、插值：在这里，就是对 $(x-\alpha_i)^r$ 建乘积树和余式树，时间为 $\widetilde O(m+D)$。([arXiv][3])

为了把实现说明完整，再展开一下恢复过程。设 $H_i=(x-\alpha_i)^K$，$H=\prod_iH_i$。计算 $z_i=(H/H_i)^{-1}\bmod H_i$ 后，答案就是 $\sum_i(r_i z_i\bmod H_i)(H/H_i)$。这些项可以沿乘积树合并。求 $z_i$ 时，把 $\alpha_i$ 平移到 $0$，就变成普通形式幂级数求逆，再平移回来即可。

矩阵 $A$ 也不必逐项展开 $R_d$。对于每个 $d$，计算所有 $S_d(k)$ 等价于求一组形如 $\sum_b c_b/(k+b+1)$ 的值。将 $c$ 倒序后与逆元序列卷积，就能一次得到全部结果，再乘上 $\frac{(k+d)!}{k!}$。所有对角线合计需要 $\widetilde O(D^2)$。

最后撤销阶乘换基。对于 $V$ 部分，先用一次卷积求出 $h_b=\int_0^1y^bf(y),dy=\sum_kf_k/(k+b+1)$，再将 $\sum_bV_{a,b}h_b$ 加到答案的第 $a$ 项，代价为 $\widetilde O(m+D^2)$。

所以，把一个重量为 $D$ 的块应用到一个长度为 $m$ 的多项式上，只需 $\widetilde O(D^\omega+mD^{\omega-2})$ 时间。 相比上一版，长多项式部分从 $mD^{(\omega-1)/2}$ 降成了 $mD^{\omega-2}$。

块的构造、总复杂度与实现

还需要在 $\widetilde O(D^\omega)$ 时间内构造块的 $U,V$，不能在这里退回立方算法。

对两个重量分别为 $D_1,D_2$ 的块，令 $D=D_1+D_2$。分别建立它们在单项式基下的截断矩阵，矩阵阶数取不小于 $2D$ 的二次幂，然后做矩阵乘法得到复合结果。

只使用结果的前 $D$ 列。这样的截断是安全的：输入次数小于 $D$，经过第一个块后次数仍然小于 $2D$，因此不会丢掉以后可能通过边界项回到低次的贡献。

从矩阵恢复核只需要 $\widetilde O(D^2)$。

先恢复 $U$。对每个 $d=1,\ldots,D$，取 $k=D-d,\ldots,D-1$。由于输出次数 $k+d\ge D$，不会受到 $V$ 的影响，所以可以直接读出 $R_d(k)=M_{k+d,k}\frac{(k+d)!}{k!}$。

这是次数小于 $d$ 的多项式在 $d$ 个连续点上的值。通过连续点插值，得到它在 $-d,\ldots,-1$ 上的值，再利用 $R_d(-b-1)=(-1)^b b!(d-1-b)!,U_{d-1-b,b}$ 恢复 $U$。连续点之间的求值转换可以由拉格朗日公式加卷积完成。

减去 $U$ 的贡献后，矩阵左上角的 $D\times D$ 部分满足 $W=VH$，其中 $H_{i,j}=1/(i+j+1)$。这个 Hilbert 矩阵满足 $H^{-1}=\operatorname{diag}(c)H\operatorname{diag}(c)$，其中 $c_i=(-1)^i(D+i)!/((D-i-1)!(i!)^2)$。因此逐行做逆元序列卷积，就能恢复 $V$。

这样，合并两个块的代价为 $\widetilde O(D^\omega)$。对操作序列分治合并，整个块的构造代价仍为 $\widetilde O(D^\omega)$。

现在取分块阈值 $B$。沿每条重链贪心划分，总重量不超过 $B$ 的操作组成普通块；单个重量超过 $B$ 的操作单独处理；每条链最后剩下的短段直接递推。

除最后的短段外，需要处理的块数为 $O(n\log n/B)$。这是因为，相邻普通块的重量之和大于 $B$，可以配对收费；被大操作隔断的短块则向该大操作收费。

所有普通块的构造和应用总代价为 $\widetilde O(nB^{\omega-1}+n^2B^{\omega-3})$。单个大操作用 NTT 乘法和一次积分完成，总代价为 $\widetilde O(n^2/B)$。

每条链最后的短段至多包含 $B$ 次操作。所有重链链头的子树大小之和为 $O(n\log n)$，因此这部分总代价为 $\widetilde O(nB)$。轻儿子的多项式乘积按大小合并，总代价为 $\widetilde O(n)$。

综合起来，

$$ T(n)=\widetilde O\!\left( nB^{\omega-1}+n^2B^{\omega-3}+nB+\frac{n^2}{B} \right). $$

取 $B=\lceil\sqrt n\rceil$，前两项同时变成 $\widetilde O(n^{(\omega+1)/2})$，后两项只有 $\widetilde O(n^{3/2})$，于是得到开头的复杂度。

块矩阵占用 $O(B^2)=O(n)$ 空间；Hermite 乘积树占用 $O(n\log n)$ 空间，因此总空间为 $O(n\log n)$。

下面摘出本次改进的核心函数。它先计算所有局部余式，然后批量乘同一个算子矩阵，最后用 Hermite 插值恢复结果；所调用的全部函数均已包含在完整源码中。

// Kernel represents:
// T[f](x) = integral_0^x U(x,y) f(y) dy
//         + integral_0^1 V(x,y) f(y) dy.
//
// This function applies only the U part.
// The caller adds the V part separately.
Poly applyVolterra(const Kernel &kernel, const Poly &input) {
    int D = kernel.d;
    int m = int(input.size());

    int K = 1;
    while (K < D) K <<= 1;

    int outSize = m + D;
    int points = 1;
    while (points * K < outSize) points <<= 1;

    // Factorial change of basis.
    Poly borel(m);
    for (int i = 0; i < m; ++i)
        borel[i] = mul(input[i], fac[i]);

    // The common matrix L:
    // degree < 2K -> degree < 3K.
    // Store it as six K-by-K blocks.
    array<Matrix, 6> A;
    for (Matrix &a : A)
        a.assign(K * K, 0);

    for (int d = 1; d <= D; ++d) {
        Poly diagonal(d);

        for (int b = 0; b < d; ++b)
            diagonal[b] = kernel.U[(d - 1 - b) * D + b];

        // v[k] = sum_b diagonal[b] / (k+b+1).
        Poly v = moments(diagonal, 2 * K);

        for (int k = 0; k < 2 * K; ++k) {
            int row = k + d;
            int blockId = (row / K) * 2 + k / K;
            int position = (row % K) * K + k % K;

            A[blockId][position] =
                mul(v[k], mul(fac[row], ifac[k]));
        }
    }

    // Each input is a remainder in the GLOBAL monomial basis:
    // borel(x) mod (x-center)^(2K).
    vector<Poly> inputs;
    {
        HermiteTree tree(points, 2 * K);
        inputs = tree.evaluate(borel);
    }

    vector<Poly> outputs(points);

    for (int first = 0; first < points; first += K) {
        array<Matrix, 2> E;
        for (Matrix &e : E)
            e.assign(K * K, 0);

        int take = min(K, points - first);

        for (int c = 0; c < take; ++c) {
            const Poly &p = inputs[first + c];

            for (int k = 0; k < int(p.size()); ++k)
                E[k / K][(k % K) * K + c] = p[k];
        }

        vector<Poly> local(take, Poly(3 * K));

        // Apply the common matrix to K local polynomials at once.
        for (int rowBlock = 0; rowBlock < 3; ++rowBlock) {
            Matrix x =
                matrixMultiply(A[rowBlock * 2], E[0], K);
            Matrix y =
                matrixMultiply(A[rowBlock * 2 + 1], E[1], K);

            for (int i = 0; i < K; ++i)
                for (int c = 0; c < take; ++c)
                    local[c][rowBlock * K + i] =
                        add(x[i * K + c], y[i * K + c]);
        }

        // Reduce L[p] modulo (x-center)^K.
        for (int c = 0; c < take; ++c) {
            int center = first + c;

            Poly p = taylorShift(local[c], center);
            p.resize(K);

            outputs[center] =
                taylorShift(p, sub(0, center));
        }
    }

    inputs.clear();
    inputs.shrink_to_fit();

    // Reconstruct L[borel(input)] from the local remainders.
    HermiteTree tree(points, K);
    Poly result = tree.interpolate(outputs);
    result.resize(outSize);

    // Undo the factorial change of basis.
    for (int i = 0; i < outSize; ++i)
        result[i] = mul(result[i], ifac[i]);

    return result;
}