多项式牛顿迭代

【前置知识】

• 泰勒展开。设 \(g\) 是一个光滑的函数，\(g(y)=\sum_{n\ge 0} \frac{g^{n}(y_0)}{n!}(y-y_0)^n\).

• 多项式 exp。

  给定多项式 \(a(x)\) 满足 \(a_0=0\)，求 \(\exp a(x)\bmod x^n\)。

  设 \(\exp a(x)=f(x)\)，设 \(f(x)\bmod x^n=f_0(x)\) 已知，试求 \(f(x)\bmod x^{2n}\)。

  根据泰勒展开公式，\(g(f(x))=g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))+\sum_{n\ge 2}\frac{f^{n}(x)}{n!}(f(x)-f(x_0))^n\)

  因为 \(f(x)\bmod x^n=f(x_0)\)，所以 \((f(x)-f(x_0))\bmod x^n=0\)，所以当 \(n\ge 2\) 时 \((f(x)-f(x_0))^n\bmod x^{2n}\) 必为 \(0\)。因此在模 \(x^{2n}\) 意义下后面的求和式可以忽略。

  于是 \(g(f(x))=g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))\pmod {x^{2n}}\)。

  取 \(g(x)=\ln x\)，上式改写为 \(\ln f(x)=\ln f_0(x)+\frac{1}{f_0(x)}(f(x)-f_0(x))\)。

  因 \(f(x)=\exp a(x)\)，则 \(\ln f(x)=a(x)=\ln f_0(x)+f_0^{-1}(x)(f(x)-f_0(x))\)。

  移项之后得到 \(f=(a-\ln f_0+1)f_0\pmod{x^{2n}}\)，可以倍增求。

【牛顿迭代】

【原版】

给定 \(g(x)\) 可导，给定 \(a(x)\)，找 \(f(x)\) 使得 \(g(f(x))=a(x)\pmod{x^n}\)。

我们可以延续上面的倍增思路，类似得到：

\[a(x)=g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))\pmod {x^{2n}}\\ f(x)=f_0(x)-\dfrac{g(f_0(x))-a(x)}{g'(f_0(x))}\pmod {x^{2n}} \]

如果 \(g'\) 好算，就可以直接倍增了。

这个一般化的模型有很多应用，先看几个基础的：

1. 取 \(g(x)=\ln(x)\)，得到 \(f(x)=\exp a(x)\)。

2. 取 \(g(x)=x^2\)，得到 \(f=a^{\frac{1}{2}}\)。

3. 取 \(g(x)=x^{-1}\)，得到 \(f=a^{-1}\)。

   代入式子看看。

   \[\begin{aligned}f(x)&=f_0(x)-(\dfrac{g(f_0(x))-a(x)}{-f_0(x)^{-2}})\\ &=f_0(x)+(g(f_0(x))-a(x))\cdot f_0^2(x)\\&=f_0(x)+f_0^2(x)(\dfrac{1}{f_0(x)}-a(x))\\&=f_0(x)+f_0(x)-a(x)f_0^2(x)\\&=f_0(x)(2-a(x)f_0(x))\pmod {x^{2n}}\end{aligned} \]

   我们发现这就是之前倍增求逆的式子。

---

这种算法为什么叫牛顿迭代？因为它和 Newton's method 的思想很像。

那 Newton's method 是什么？我们看个例子。

求 \(x^2=w\) 的 \(x\)。（求 \(\sqrt w\)）

初始随便猜一个 \(x_0\)。然后 \(x_{k+1}=N(x_k)=x_k-\dfrac{(x_k^2-2)}{2x_k}\)。

即令 \(\phi(x)=x^2-2\)，\(x_{k+1}=x_k-\dfrac{\phi(x_k)}{\phi'(x_k)}\)。迭代无穷次。

这是为什么呢？从函数图像上考虑，就是取 \(x_k\) 处斜率直线与 \(x\) 轴交点作为 \(x_{k+1}\).

---

回顾上面多项式倍增法，其实是类似的思路。是定义了 \(\phi(f(x))=g(f(x))-a(x)\pmod {x^n}\)。

这里怎么理解？\(\phi\) 是为 \(k[x]\rightarrow k[x]\) 的映射，\(k[x]\) 是所有实系数多项式。问题就是找 \(\phi(f(x))=0\) 的根 \(f(x)\)。（注意不是根 \(x\)！！）\(a(x)\) 是常值多项式，也就是我们不将其视作多项式，而视作一个常数。（类比求 \(kx+b=0\) 的根）

回顾上面的公式：

\[f(x)=f_0(x)-\dfrac{g(f_0(x))-a(x)}{g'(f_0(x))}\pmod {x^{2n}} \]

发现 \(g(f_0(x))-a(x)=\phi(f_0(x))\)，\(g'(f_0(x))=\phi'(f_0(x))\)，注意 \(a(x)\) 是常数，求导之后没有了。因此可以改写为 \(f(x)=f_0(x)-\frac{\phi(f_0(x))}{\phi'(f_0(x))}\)。

所谓的 "多项式倍增法" 就是原版牛顿迭代从实数拓展到多项式了。

但我们还能拓展！！！

【拓展】

原问题：求 \(\phi(f(x))=g(f(x))-a(x)\) 的根 \(f(x)\) 使得 \(\phi(f(x))=0\)。

因为 \(\phi(f(x))=g(f(x))-a(x)\)，所以就是求 \(g(f(x))=a(x)\).

ex（在线版）：\(a(x)\) 初始未给定全部系数，而是与 \(f(x)\) 系数相关（类似自卷）。

【例一】

（为了方便，直接取 \(g=\ln\) 了）

给定 \(g,c\)。要求 \(f(x)=(1-a(x))\cdot c(x)\)，保证 \(c_0=1\)，\(a_0=0\)，\(f_0=1\)。

求 \(\ln(f(x))=a(x)\) 的根 \(f(x)\)。

---

换句话说，就是要找两个函数 \(a,f\)，既满足 \(f=(1-a)c\)，又满足 \(\ln f=a\)。

再换句话说，要找一个函数 \(a\)，满足 \(\ln((1-a)c)=a\iff (1-a)c=\exp a\)。

再再换句话说，要找一个函数 \(a\)，满足 \(\dfrac{\exp a}{1-a}=c\)。

原本的形式只能解决 \(\exp a=c\)，这里我们拓展了一个 \(1-a\) 分母。

那么这个问题怎么解呢？

首先定义 \(\phi(f(x))=\ln (f)-a\)。根据牛顿迭代：

\[f(x)=f_0(x)-\dfrac{\phi(f_0(x))}{\phi'(f_0(x))}\pmod{x^{2n}} \]

所以 \(f=f_0(1+a-\ln f_0)\)。代入 \(f=(1-a)c\)：\((1-a)c=f_0(1+a-\ln f_0)\)。

移一下，\((c-f_0+f_0\ln f_0)=a(c+f_0)\pmod{x^{2n}}\)。因为 \(c_0=f_0=1\)，所以 \(c+f_0\) 可逆，所以 \(a=(c-f_0+f_0\ln f_0)(c+f_0)^{-1}\pmod {x^{2n}}\)。右侧在倍增时全部已知，可以求 \(a\)，然后 \(f=(1-a)c\) 也可以求出来了。

复杂度 \(T(n)=T(n/2)+n\log n=n\log n\)。

是不是感觉很对？上面的做法是错的！！！

虽然式子正确，但是 \(a(x)\bmod x^{2n}\) 未知！（\(a(x)\bmod x^n\) 已知）

那怎么做？

---

这种在线的卷也有 cdq 分治解法，两个 \(\log\)，但是因为常数，不一定慢。

牛顿迭代的主要优势在于推导简单。

【例二】

给定 \(c\)，求 \(a,f\)：\(f=\dfrac{c}{1+a}\) 且 \(\ln f=a\)。

等价一下。

已知 \(c\)，求 \((1+a)\cdot \exp a=c\)。（保证 \(a_0=0,c_0=1\)）

定义 \(\phi\) 是一个自变量多项式的函数，\(\phi(f(x))=\ln f(x)-a(x)\).

类似地，\(f=f_0(1-\ln f_0+a)\pmod {x^{2n}}\)，代入 \(f=\dfrac{c}{1+a}\)。

代入得到 \(c=(1+a)(f_0(1+a-\ln f_0))\)。

换元，\(b=1+a\)，显然求 \(a\) 和求 \(b\) 是等价的。\(c=bf_0(b-\ln f_0)\)。

这个方程可以写成 \(b^2-p\cdot b-q=0\) 的形式，且 \(p_0=0\)。（注：\(f_0\) 有逆）

采用最朴素的配方法，\((b-\frac{p}{2})^2=q+(\frac{p}{2})^2=d\)。\(p,q,d\) 都是已知的。

这是有唯一解的。因为 \(b_0=a_0+1=1\)，\(d\) 已知，可以一步一步推出来。

所以用一个多项式开方可以 \(O(n\log n)\) 求 \(b\)，\(a=b-1\)，进而再求出 \(f\)。

所以一次迭代还是 \(O(n\log n)\) 的。\(T(n)=T(n/2)+O(n\log n)=O(n\log n)\)。

【例三】

来点实际问题。

loj 6538 烷基计数 加强加强版

\(n\) 个结点的无标号有根树计数，要求每个结点的子结点个数 \(\le 3\)。认为子树无序。

\(n\le 10^5\)，答案 \(\bmod \;998244353\)。

记答案为 \(f_n\)。首要肯定是找 \(f\) 的递推式。

一个 naive 的想法是 \(f_n=\sum_{a+b+c=n-1}f_af_bf_c\)，这显然是会算多的，因为子树有序了。

另一个 naive 的想法是 \(f_n=\sum_{a+b+c=n-1}\frac{f_af_bf_c}{6}\)。这是会算少的，因为只有 \(a,b,c\) 互不相等才除以 \(6\)，其他情况不用除以 \(6\) 这么多。

那正确的想法是什么呢？考虑容斥。

令 \((X,Y,Z)\) 有序三元组表示一颗子树有序的树：根的三颗子树按顺序分别是 \(X,Y,Z\)。

注：要求 \(|X|+|Y|+|Z|=n-1\).

然后呢？枚举 \(X,Y,Z\) 这三棵树的形态太恶心了，我们只需要枚举 \(|X|,|Y|,|Z|\) 即可。

先让 \(f_n\leftarrow\sum_{a+b+c=n-1}f_af_bf_c\cdot \frac{1}{6}\)。此时第一类算的正好，第二类和第三类不对。

考虑第二类，它形如 \((X,X,Y),(X,Y,X),(Y,X,X)\)。这三种在上面会各算一次，也就是这种方案一共算了 \(3\times \frac{1}{6}=\frac{1}{2}\) 次，也就是 \(X,X,Y\) 的排列已经贡献了 \(\frac{1}{2}\)。但是我们希望 \(X,X,Y\) 的排列一共贡献 \(1\)，差了 \(\frac{1}{2}\)。所以我们加上一点东西：

\(f_n\leftarrow \sum_{2a+b=n-1}f_af_b\cdot \frac{1}{2}\)。加上这一部分之后，第一类第二类都算好了。只剩第三类了。

对于 \(X,X,X\)，它在 \(f_af_bf_c\) 的时候会统计 \(1\) 次，贡献 \(\frac{1}{6}\)；在 \(f_af_b\) 的时候会统计 \(1\) 次，贡献 \(\frac{1}{2}\)。也就是它在先前的统计里已经做了 \(\frac{2}{3}\) 的贡献，还差 \(\frac{1}{3}\)。

因此 \(f_n\leftarrow \sum_{3a=n-1}f_a\cdot \frac{1}{3}\)。

综合起来，\(f_n=\sum_{a+b+c=n-1}\frac{1}{6}f_af_bf_c+\sum_{2a+b=n-1}\frac{3}{6}f_af_b+\sum_{3a=n-1}\frac{2}{6}f_a\)。

于是得到一个 \(O(n^4)\) 的暴力方法，然而 \(n\le 10^4\).

因为右边有 \(a+b+c=n-1\) 这种形式的东西，考虑母函数优化。

\[f(x)=1+x\cdot \dfrac{f(x)^3+3f(x^2)f(x)+2f(x^3)}{6} \]

然后思考一下怎么解，直接上倍增。

设 \(f_0=f\bmod x^n\) 已知，目标求 \(f\bmod x^{2n}\)。

\(f_0\) 已知，则 \(f(x^2),f(x^3)\bmod x^{2n}\) 的系数已知。那么我们记 \(B(x)=f(x^2)\bmod x^{2n},C(x)=f(x^3)\bmod x^{2n}\)。

改写方程：

\[f=1+x\cdot \dfrac{f^3+3Bf+2C}{6}\pmod {x^{2n}} \]

省略一些移项的 dirty work，得到：

\[x\cdot f^3+(3xB-6)\cdot f+(2xC+6)=0\pmod {x^{2n}} \]

可以写作 \(xf^3-pf+q=0\) 的形式。考虑牛顿迭代。

定义 \(\phi(f)=x\cdot f^3-p\cdot f+q\)。有 \(f=f_0-\frac{\phi(f_0)}{\phi'(f_0)}\pmod {x^{2n}}\)。

\(f_0\) 是已知的。要解 \(f\)。

\(\phi(f_0)\) 直接代入即可。\(\phi'\) 呢？\(\phi'(f)=3x\cdot f^2-p\)。取 \(f_0\) 点值可以算 \(\phi'(f_0)\)。

【总结】

牛顿迭代的重点是要知道 \(\phi'\) 是把 \(f\) 整个多项式当作自变量求导。