牛顿迭代实现多项式半家桶及其常数优化——2025.6.15鲜花

牛顿迭代实现多项式半家桶及其常数优化——2025.6.15鲜花

ネバーランド

なんでそんなこと言い出すの
あたし何も悪くないよ
きみが好きを始めたのに
あたしそれを信じたのに
きみのためになんでもしたのに
きみ以外で幸せ
なれっこないんだもん
そこんところ
きみはどう思ってるの
さあ一緒にいこう？ネバーランド
何度でもきみと遊びたい
バイバイ言ったら死刑なルール
絶対誰にも渡さない
嗚呼こんなに甘いネバーランド
どうしてきみはつらそうなの？
どんまい？ 飽いたらリセットルール？
どうしてもきみが離れない
まだ好きだよね？
バッドなエンドだ あっはっは
無理すぎて
あれもこれも病んじゃった
キュートな愛も無駄
「重すぎてしんどい」って
泣いちゃった
バッドなエンドだ あっはっは
無理すぎて
あれもこれも病んじゃった
キュートな愛も無駄
「重すぎてしんどい」って
泣いちゃった
変な冗談は良くないかも
嘘と言ってお互いのために
きみが好きを始めたのに
あたしそれを信じたのに
きみのためになんでもしたのに
あたし以外で幸せ
なれっこないでしょ？
ありえないよ 笑えない
許してやんない
さあ一緒にいこう？ネバーランド
「最後」とかなんで？認めない
頂戴でっかいごめんのポーズ
なんとかやり直せるんじゃない？
終わっちゃうくらいなら
嫌いにさせてよ
これ以上
あたしに優しくしないでよ
だって好きって思ってしまうもん
ずっと一緒願ってしまうもん
やっぱ全部きみが悪いんだ
じゃあ一緒に死のう？ネバーランド
何度でもきみで拗れたい
最強ふたりの未来のルーム
絶対誰にも渡さない
嗚呼こんなに甘いネバーランド
やっときみも笑えるんだね
天才 あたしの想いがルール
最後に愛は勝つ そうでしょう？
ハッピーエンドだ あっはっは
これからもきみがすきなんだった
嫉妬も酸素もバイバイ
これからもだいすきなんだった
ハッピーエンドだ あっはっは
これからもきみがすきなんだった
嫉妬も酸素もバイバイ
これからもだいすきなんだった
あは

卡了一天半的常的结果，其实挺好玩的。

首先会一下牛顿迭代。

我们有一个非平凡函数 \(G(x, y)\)，求模 \(x ^ n\) 意义下的 \(f(x)\) 满足 \(G(x, f(x)) \equiv 0 \pmod {x ^ n}\)

我们先把 \(n\) 补齐到 \(2 ^ k\)。

考虑倍增，若我们已经求出模 \(x ^ n\) 意义下的 \(f_0\)，考虑求出模 \(x ^ {2n}\) 意义下的 \(f\)。

将 \(G(x, f)\) 在 \(f = f_0\) 处泰勒展开，有：

\[\sum_i \frac{\frac{\mathrm d^i G}{\mathrm d^i y} (x, f_0)}{i!}(f - f_0) ^ i \equiv 0 \pmod{x ^ {2n}} \]

考虑到 \(f - f_0\) 的最低次项是 \(n\)，则 \(i \ge 2\) 时其值模 \(x ^ {2n}\) 为 \(0\)。所以原式等价于：

\[G(x, f_0) + \frac{\mathrm d G}{\mathrm d y}(x, f_0)(f - f_0) \equiv 0 \pmod{x ^ {2n}} \]

\[f \equiv f_0 - \frac{G(x, f_0)}{\frac{\mathrm d G}{\mathrm d y}(x, f_0)} \pmod{x ^ {2n}} \]

然后就做完了。

具体的式子可以看 oi-wiki

考虑优化常数，这里大量参考 negiizhao 的 uoj 博客，但只写了其中比较浅显的一点，更深的可以看原文章。

以下统一设 \(k_0\) 表示上次迭代 \(k\) 的结果，\(k\) 表示本次迭代或真实值，\(\hat{g^n}\) 表示 \(g\) DFT 以后的长为 \(n\) 的点值，将 DFT 操作和 IDFT 操作统称为 FFT 操作，简称 FFT。

实际实现中存在很多的 FFT 复用，这里我们假设一个 FFT 的结果求出来后就会一直使用。

设一次大小为 \(n\) 的 FFT 的代价是 \(E\)，将 FFT 的常数视为线性的。

由于我的疏忽，可能会有算错的，还请指正。

1. 求逆，\(g = \frac 1f\)：

   式子：\(g = 2g_0 - g_0^2f\)

   容易发现对于 \(g_0^2f\) 我们要做一次大小是 \(2n\) 的卷积，一次大小是 \(4n\) 的卷积，总计 \(18E\)。

   我们有技巧可以用 \(2E\) 的代价从 \(\hat{g^n}\) 求 \(\hat{g^{2n}}\)，从而使代价降到 \(16E\)。这个技巧将会在 Bostan-Mori 的时候讲到，这里我们有更好的办法。

   发现 \(g_0f\) 的前 \(n\) 项只有常数项是 \(1\)，所以我们拆一下式子：\(g = 2g_0 - g_0^2f = g_0 - (fg_0 - 1)g_0\)，这样 \(fg_0 - 1\) 的系数在 \([n, 3n)\)，则可以直接做 \(2n\) 的循环卷积，同样的 \((fg_0 - 1)g_0\) 也可以做 \(2n\) 的循环卷积，注意中间有清空操作，所以总代价为 \(10E\)。

2. 开根，\(g = \sqrt f\)：

   式子：\(g = \frac {g_0}2 + \frac f{2g_0}\)，暴力做需要 \(2n\) 的求逆和 \(4n\) 的卷积，总计 \(22E\)。

   类似求逆的优化，我们吧式子展开成 \(g_0 - \frac{g_0 ^ 2 - f}{2g_0}\)，此时 \(g_0 ^ 2 - f\) 是 \(2n\) 的循环卷积，\(\frac{g_0 ^ 2 - f}{2g_0}\) 依然是 \(2n\) 的循环卷积。注意到 \(g_0 ^ 2 - f\) 最低次项是 \(n\)，则 \(\frac 1{2g_0}\) 只需要求 \(n\) 项，总计 \(18E\)。

   我们接着发现，\(g_0 ^ 2 - f\) 的系数竟然是 \([n, 2n)\)，也就是说可以换成长为 \(n\) 的循环卷积，总计 \(13E\)。

   注意这个结果要将求逆展开来和开根一起迭代进行 FFT 结果的复用。

   注意到开根迭代的 FFT 结果中有一个可以和下一轮迭代复用，最终是 \(12E\)。

3. 求商或求 \(\ln\)，\(h = \frac fg\)：

   \(\ln f = \int \frac{f'}f\)，其中求导和积分都是 \(\mathcal{O}(n)\) 的，所以和求商等价。

   暴力做是求逆加 \(4n\) 的卷积，总计 \(22E\)

   依然是展开：\(h = h_0 - \frac{gh_0 - f} f\)，类似开根用循环卷积优化，代价是 \(15E\)，类似的，需要将求逆展开来和开根一起迭代。

   有一个技巧是通过 \(\hat{g_0}\) 快速求出 \(fg\)，具体就是将 \(g\) 拆开，懒得写了。

4. 指数，\(g = \exp f\)

   式子：\(g = g_0(1 - \ln g_0 + f)\)。

   暴力做是 \(\ln\) 加 \(4n\) 卷积，总计 \(27E\)

   考虑拆式子，我们拆成 \(g = g_0 - g_0\int{(\frac{g_0'}{g_0}-f')}\)

   我们发现这次的 \(\frac 1{g_0}\) 是在 \(2n\) 次截断，所以没法将求逆展开来和开根一起迭代，只能暴力做，但是循环卷积依然一样做，代价是 \(22E\)。

   我们继续推式子，设 \(h_0 = \frac 1{g_0} \bmod x ^ n\)：

   \[\begin{aligned} g \bmod x^{2n} &= g_0 - g_0(\int{(\frac{g_0'}{g_0})}-f) \bmod x^{2n} \\ &= g_0 - g_0\int{(\frac{g_0'}{g_0}-f')} \bmod x^{2n} \\ &= g_0 - g_0\int{(g_0h_0(\frac{g_0'}{g_0}-f'))} \bmod x^{2n} \\ &= g_0 - g_0\int{(g_0'h_0-f'-(g_0h_0-1)f_0')} \bmod x^{2n} \end{aligned}\]

   这里面较为关键的一步是发现 \(\frac{g_0'}{g_0}-f'\) 的系数 \(\ge n - 1\)，则 \(g_0h_0(\frac{g_0'}{g_0}-f') \equiv \frac{g_0'}{g_0}-f' \pmod {x ^ {n - 1}}\)，积一下就是 \(\int(g_0h_0(\frac{g_0'}{g_0}-f')) \equiv \int(\frac{g_0'}{g_0}-f') \pmod {x ^ n}\)。

   这里的 \(g_0'h_0-f'\) 和 \(g_0h_0-1\) 都可以 \(n\) 循环卷积，\((g_0h_0-1)f_0'\) 是 \(n ^ 2\) 循环卷积 \(g_0\int{(g_0'h_0-f'-(g_0h_0-1)f_0')}\) 是 \(n ^ 2\) 循环卷积，将求逆展开来和开根一起迭代，依然有一个 FFT 结果可以对下轮复用。

   发现 \(g_0h_0-1 \to (g_0h_0-1)f_0'\) 这一步中有一个将前一半多项式移到后一半，这个可以先二倍点值在暴力吧点值乘上 \(x^k\)，代价可以做到 \(2E\)。

   如何二倍点值见 Bostan-Mori 部分。

   实现下来代价是 \(19E\)

5. Bostan-Mori

   本段参考 EI 一种新的线性递推计算方法

   先会一下 Bostan-Mori，求 \([x^k]\frac {F(x)}{Q(x)}\)，考虑：

   \[[x^k]\frac {F(x)}{Q(x)} = [x^k]\frac {F(x)Q(-x)}{Q(x)Q(-x)} \]

   我们发现 \(Q(x)Q(-x)\) 只有偶次项，而 \(F(x)Q(-x)\) 随 \(x\) 的奇偶只会有奇次项或偶次项有用，于是折半即可，设多项式长度是 \(n\)，容易发现多项式长度不会变。复杂度是 \(n\log n\log k\)。

   朴素做代价是 \(10E\)，这里记 \(F = F(x), A = Q(x), B = Q(-x)\)

   首先的一个技巧是 \(\hat{A_i} = \hat{B_{i \oplus \frac n2}}\)，所以可以省掉 \(B\) 的一遍 \(2n\) 的 FFT。

   然后考虑如何在点值上提取奇偶次幂。以偶次为例，奇次是类似的。我们设 \(C_i = \frac {A_i + (-1) ^ i A_i}2\)，则 \(\hat{C_i} = \frac {\hat{A_i} + \hat{A_{i \oplus \frac n2}}}2\)，只要取 \(\hat{C_i}\) 的前一半即可。

   最后就是如何直接推出二倍点值，这也是容易的

   \[\hat{A^{2n}_2k} = \hat{A_k} \]
   \[\hat{A^{2n}_{2k + 1}} = \sum_i\hat{\omega_{2n}^i A_i}\omega_{2n}^{ik} \]

   因此我们先 IDFT 得到 \(A\)，再乘 \(\omega_{2n}^{i}\)，在做 DFT 即可。

   最后就只有 \(4E\)。

具体的实现可以看我板子（写的很屎）。

P