草稿。。。

终于做出来一点原创的结果了。

上一篇写了 4 阶牛顿迭代，最后提到还有一个 18E 的 pow 想法。现在也写出来。

问题定义

多项式幂的定义是

$$g = f^k \bmod {x^n}, \qquad f(0)=1$$

最直接的做法是

$$g = \exp(k \ln f).$$

一般来说，如果你使用的是比较朴素的 13E ln 和 17E exp，得到的就是一个 30E 的 pow 算法；如果你选择的是 SOTA 10E ln + 11.5E exp，你将收获一个 21.5E 的 pow；根据我的实践，10E ln + 14E exp 在实机上表现最好，也就是 24E 的 pow。

也有独立的 pow 算法，目前的 SOTA 是 20.25E，本文提出了一种新的 18E 的计算多项式幂的算法。

之前尝试在线卷积的时候，就有想过能不能同时维护两边的迭代式，得到更好的复杂度，本次就继续沿着这个思路取得了一些进展。

分块

建议阅读 多项式牛顿迭代的分块优化，本文引用了其不少结论。

记块长为 $m$，$X=x^m$，把

$$f = f_{[0]} + f_{[1]}X + \cdots,\qquad g = g_{[0]} + g_{[1]}X + \cdots.$$

假设已经知道 $g_{[0]},\cdots,g_{[s-1]}$，现在要求 $g_{[s]}$。

在前一篇 inv 的分块中，我们有一个结论，首先记 $fg$ 在缺少 $g_{[s]}$ 下的卷积为 $\psi$，即

$$\psi = (f_{[0..s]}g_{[0..s-1]})_{[s]}$$

我们有一个等式：

$$\begin{aligned} (fg)_{[s]} - \psi &= (f_{[0..s]}g_{[0..s]})_{[s]} - (f_{[0..s]}g_{[0..s-1]})_{[s]} \\ &= (f_{[0..s]} g_{[s]} X^{s})_{[s]} \\ &= f_{[0]} \ast g_{[s]} \bmod X \end{aligned} \tag{1}$$

接下来我们将在两边的迭代式考虑这两个卷积的差。

微分方程

沿用之前的记号，考虑原地微分算子 $\Delta h = xh'(x)$。

为了方便起见，我们记 $\Delta_i f = X^{-i} \Delta(X^i f)$，也就是一个平移到微分算子，有

$$\Delta (f_{[0]} + f_{[1]} X + \cdots) = (\Delta_0 f_{[0]}) + (\Delta_1 f_{[1]}) + \cdots$$

考虑 $g=f^k$，两边微分得到

$$f\Delta g = k g \Delta f$$

把两边同时加上一个 $g \Delta f$，写成

$$\Delta(fg) = (k+1)g\Delta f$$

这个式子更适合分块，因为左边出现的是 $fg$，而最后我们要解的是 $g$ 的新块。其分块形式为

$$(\Delta fg)_{[s]} = (k+1)(g \Delta f)_{[s]} \tag{2}$$

考察两侧式子缺少第 $s$ 块的情况，我们考虑维护左侧 $fg$，考虑应用公式 (1)，得到

$$(fg)_{[s]} - \psi = f_{[0]}g_{[s]} \bmod X.$$

对于右侧，考虑维护 $g \Delta f$，继续应用公式 (1)，得到

$$(g \Delta f)_{[s]} - \phi = (\Delta f_{[0]}) g_{[s]} \bmod X$$

再考虑等式 (2)，因此有

$$\Delta_s (\psi + f_{[0]}g_{[s]}) = (k+1)(\phi + g_{[s]}\Delta f_{[0]}) \pmod X$$

接下来要把 $g_{[s]}$ 解出来。直接求解比较困难，两侧一边是 $g_{[s]}$ 另一边是 $\Delta_s g_{[s]}$，消不掉。

考虑换元，令

$$g_{[s]} = g_{[0]}u.$$

因为 $g_{[0]} = f_{[0]}^k$，所以

$$\Delta (f_{[0]}g_{[0]}) = f_{[0]}\Delta g_{[0]} + g_{[0]}\Delta f_{[0]} = (k+1)g_{[0]}\Delta f_{[0]}.$$

再注意到

$$\begin{aligned} \Delta_s(f_{[0]}g_{[s]}) &= \Delta_s(f_{[0]}g_{[0]}u) \\ &= f_{[0]}g_{[0]}\Delta_s u + u\Delta (f_{[0]}g_{[0]}) \\ &= f_{[0]}g_{[0]}\Delta_s u + (k+1)g_{[0]} u \Delta f_{[0]} \\ &= f_{[0]}g_{[0]}\Delta_s u + (k+1)g_{[s]} \Delta f_{[0]} \end{aligned}$$

这恰好是公式 (2) 中的两项，因此

$$f_{[0]}g_{[0]}\Delta_s u = (k+1)\phi - \Delta_s\psi.$$

也就是

$$u = \Delta_s^{-1}\left(\frac{(k+1)\phi - \Delta_s\psi}{f_{[0]}g_{[0]}}\right), \qquad g_{[s]} = g_{[0]}u.$$

具体过程

计算 $\psi$ 和 $\phi$ 的方法与上一篇分块的方法完全相同，即计算 $fg$ 和 $g \Delta f$ 的第 $s$ 块。由于 $g$ 是共享的，只需要维护 $f,\Delta f,g$ 的分块 DFT，再通过两次 IDFT 即可获得 $\psi$ 和 $\phi$。

具体计算流程如下：

0. 假设已经求出求出 $g_{[0]} = f_{[0]}^k$。
1. 计算 $\mathcal{F}_{2m}(g_{[0]})$，这里 2E 的开销被计入迭代内，故无开销。
2. 计算 $f_{[0]}g_{[0]} \bmod X$，需要 4E。
3. 预处理 $h = (f_{[0]}g_{[0]})^{-1}$ 可以使用任何 6E 的算法。
4. 依次计算第 $s$ 块：
   1. 计算 $\mathcal{F}_{2m}(f_{[s]})$ 和 $\mathcal{F}_{2m}(\Delta_s f_{[s]})$ 和 $\mathcal{F}_{2m}(g_{[s-1]})$，开销 6E。
   2. 用块乘积算出 $\psi$，开销 2E。
   3. 用块乘积算出 $\phi$，开销 2E。
   4. 计算 $\rho = (k+1)\phi - \Delta_s\psi$，线性开销。
   5. 计算 $v = h\rho \bmod X$，开销 4E。
   6. 计算 $u = \Delta_s^{-1}v$，线性开销。
   7. 计算 $g_{[s]} = g_{[0]}u \bmod X$，开销 4E。

总之，单个新块的复杂度为 18E，当块数 $r \to \infty$ 时，均摊复杂度也就是 18E。

严格的说，复杂度 $T(n) \leqslant 18\mathsf{E}(n) + O(n \sqrt{\log n})$。

实现与测试

性能测试取 NTT 模数 $P=998244353$，指数固定为 $k=7$，并使用 AVX2 优化。

可以看到收益比较稳定，大概快 20% 多一点，符合预期。

总结

本来只是闲暇时的瞎想，没想到越琢磨越有道理，终有此篇。也总算是做出来一点点属于自己的东西了。