DGF 的计算理论：牛顿迭代与特殊求逆运算

飞雨烟雁 · 2024-12-11 18:31:59 · 算法·理论

人类对 DGF 的开发尚不足百分之一。

前言

本文分为两个部分。第一部分是「牛顿迭代」，在这部分我们会呈现如何将牛顿迭代运用到 DGF 之上。第二部分是「特殊求逆运算」，包括类积性函数求逆{}和广义求逆。一些记号约定及定义如下：

1. 我们用花体字母代表多项式，例如 \mathcal F,\mathcal G；用 F,G 等表示 DGF，并用对应的小写字母表示 DGF 所对应的数论函数；

2. 多元向量简写为 \mathbf x=(x_1,x_2,\cdots,x_m)；

4. 定义多元多项式 \mathcal F(\mathbf x) 的最低次数为所有项的次数的最小值。例如 xyz+x^4 的最低次数是 3；

5. 对于 DGF G=\sum_{k}a_kk^{-x}，其前 n 项指的是 a_1,a_2,\cdots,a_n，其第 n 项用 G[n] 表示。

DGF 的牛顿迭代

多元多项式的牛顿迭代

DGF 牛顿迭代的理论基础是多元多项式的牛顿迭代。

给定具有简单形式的 \mathcal F(x)，试求多元多项式 \mathcal G，使得 \mathcal F(\mathcal G(\mathbf x))\equiv 0\pmod {\mathbf x^n}。

注：此处 {}\!\bmod \mathbf x^n 的意思是去掉所有次数 \ge n 的项。

考虑倍增，假设我们已经求出了 {}\!\bmod\mathbf x^n 时的解 \mathcal G_0，我们要怎么求得 \mathcal G 使得 \mathcal F(\mathcal G(\mathbf x))\equiv 0\pmod {\mathbf x^{2n}} 呢？

思路和一元多项式的牛顿迭代是一样的，我们考虑 \mathcal F 在 \mathcal G_0 处的泰勒展开：

因为 \mathcal G(\mathbf x)-\mathcal G_0(\mathbf x) 的最低次数是 n，所以红色部分在 {}\!\bmod \mathbf x^{2n} 时为 0。因此有：

若 \mathcal F'(\mathcal G_0(\mathbf 0))\neq 0，则有：

这和一元多项式的式子一模一样。因此，若 \mathcal F\circ \mathcal G 能在 \Theta(n\log n) 时间内算出，则多元多项式的牛顿迭代可以在 \Theta(n\log n) 时间内完成。

DGF 的牛顿迭代

然后我们来看 DGF 的牛顿迭代，具体而言，我们希望解决如下问题：

设所有 DGF 构成环 R。给定具有简单形式的 \mathcal F(x)\in R[x]，求 G\in R，使得 \mathcal F\circ G=0。

先给个例子，当 \mathcal F(x)=x^2-3x+1+\zeta 时，G=\frac 32\pm\sqrt{\frac 54-\zeta}。

下面提供一个在 \Theta(n\log n\log \log n) 时间内求出 G 的前 n 项的算法。

一个自然的想法是，我们建立从 R 到 \mathbb C[x_1,x_2,\cdots,x_n,\cdots] 的双射 \rho，使得对于 G\in R，有：

注：v_p(n)=\max\{\beta:p^\beta\mid n\}。

可以发现，这是两个环之间的同构。也就是说，想要在 R 中求解 \mathcal F\circ G=0，我们只需在 \mathbb C[x_1,\cdots] 上求解 \mathcal F(\mathcal H)=0，且 G=\rho^{\left<-1\right>}(\mathcal H)。

尽管 \mathbb C[x_1,\cdots] 上有无限个未定元，但考虑到我们只关心 G[1],G[2],\cdots,G[n] 的取值，所以我们只看前 \pi(n) 个未定元。

这样，我们就把 DGF 的牛顿迭代问题转化为多元多项式的牛顿迭代问题。当然，实际计算中根本不需要 DGF 转多元多项式，我们只是借助多元多项式这个框架来说明牛顿迭代的适用性。

下面的核心议题是，我究竟要迭代多少次才能正确求出 G[1],G[2],\cdots,G[n]。根据前面的推导，每次迭代都会使正确次数翻倍，所以我们只需迭代 \Theta(\log n) 次？事实上，因为 DGF 的特殊性，迭代 \Theta(\log \log n) 就足够了。

具体地讲，想要求出 G[1],\cdots,G[n]，我们只需求出 \mathcal F(\rho(G))\equiv 0\pmod {\mathbf x^{N+1}}，其中 N=\max_{1\le i\le n}\Omega(i)，\Omega 是质因子个数（计重）。显然 \Omega 在 2 的幂次时取最大值，也即 N=\lfloor\log_2n\rfloor。那么最终只需迭代 \Theta(\log N)=\Theta(\log \log n) 次。

假设 \mathcal F\circ G 能在 \Theta(n\log n) 时间内完成计算，那么我们就得到了一个 \Theta(n\log n\log\log n) 的算法。

非齐次牛顿迭代

当然，部分读者可能会有这样的疑惑，为什么时间复杂度是 \Theta(n\log n\log \log n)，而不是 T(n)=T(\sqrt n)+\Theta(n\log n) 算出来的 \Theta(n\log n) 呢？这是因为，在上面的理论框架下，想要保证迭代的正确性，第 k 轮迭代后必须保证 \mathcal F(\rho(G))\equiv 0\pmod {\mathbf x^{2^k}}。也就是说，第一轮迭代后必须确保所有 \pi(n) 个未定元的系数是正确的。因此，即便是在第一轮迭代，我们也必须进行长度为 n 的 DGF 计算。所以这 \log N 轮迭代必须做 \log N 次长为 n 的 DGF 计算，也就是时间复杂度 \Theta(n\log n\log N)=\Theta(n\log n\log \log n)。

但是，实际计算又告诉我们，似乎在第 k 轮计算时，我们只需做长为 2^{2^k} -1 的 DGF 计算即可。为了解释这样做的正确性，我们需要新的理论框架，即非齐次牛顿迭代。

我们先用二元多项式举个例子，假设我们已知 \mathcal F(\mathcal H_0(x,y)) 中 1,x 的系数是正确的，那么牛顿迭代中的 \mathcal H-\mathcal H_0 的低次项就有 x^2,xy,y。那在平方过后，(\mathcal H-\mathcal H_0)^2 的低次项就是 x^4,x^3y,x^2 y,xy^2,y^2。换言之，按牛顿迭代的式子算出的 \mathcal H 中的 x^2,x^3,xy,y 的系数是正确的。如下图，浅绿色是已知正确部分，红色是会受到 (\mathcal H-\mathcal H_0)^2 影响的部分，深绿色是迭代出来后新增的正确部分。

根据这个想法，我们就有非齐次的牛顿迭代。假设 \mathcal H-\mathcal H_0 的低次项构成集合 A，则牛顿迭代后低于 A^2:=\{\mathbf x\mathbf y:\mathbf x,\mathbf y\in A\} 项就都是正确的。推广到 DGF 上，就是说，如果 G-G_0 的低次项构成集合 A，则低于 \{xy:x,y\in A\} 的项在迭代后都是正确的。当然，此处的低于应该理解为整除。因为多元单项式的偏序关系在同构 \rho 的意义下等价于正整数集的整除关系。

因此，假设我们求出的 G_0 能保证 G_0[1],G_0[2],\cdots,G_0[m-1] 都是正确的，那么 A=\{m,m+1,\cdots\}，即 A^2=\{m^2,m(m+1),\cdots\}，故低于 m^2 的部分都是正确的。因此我们只需取前 m^2-1 项参与 DGF 计算即可。

上面是理论介绍，实际操作很简单，这里我们给个伪代码：

DGF Newton_Iteration(int n){
    G[1] = Initial_Value;
    for(long long m = 2; m <= n; m = m * m){
        int Deg = min(m * m - 1, n);
        G = G - Composite(F, G, Deg) / Composite(Derivate(F), G, Deg);
    }
    return G;
} 

假设计算 \mathcal F,\mathcal F' 复合 G 的前 n 项的时间复杂度是 \mathsf M(n)，则该做法的时间复杂度是：

注：此处假定 \mathsf M(n)=o(n^{1+\varepsilon})。

更精确的说，最高阶的系数不会超过 2，即复杂度上界为 (2+o(1))\mathsf M(n)。

这里我们解释一下为什么系数是 2。因为当 n=2^{2^k} 时，比如 n=256，此时我们每轮计算的长度分别为 3,15,255,256。可以发现，为了计算 256 这一项，我们又花了一次时间来计算，因此系数是 2。我们当然可以进一步优化，因为计算 256 这一项只需要看 2,4,8,\cdots,128 的系数，完全可以快速计算。但这些都是后话了，无论如何，我确信我们已经达到了理论复杂度下界，后人只需要在系数上修修补补即可（

另外，进一步的优化在实践意义下是没必要的。因为一般遇到的 n 大概是 2\times 10^6 量级，而每轮迭代的长度分别为 2,4,16,256,65536,4294967296,\cdots，此时不会出现上面说的情况。

DGF 的特殊计算

类积性函数求逆（半在线卷积）

给定积性函数 g(n)，若 f(n) 满足 f(n)=\sum_{d|n,d<n}f(d)g(n/d),f(1)=1，试求出 f(1),\cdots,f(n) 的值。

换言之，我们就是要求 F=1/(2-G)，这相当于是做 DGF 求逆。回忆求逆的牛顿迭代式：

当 F_0 的前 \sqrt n 项都是正确的时，我们就能通过上面这条式子直接求出前 n 项。在上面这条式子中，计算 F_0^2 的时间复杂度是 \Theta(n)，乘上 2-G 的时间复杂度是 \Theta(n\log\log n)。至于 F_0，我们可以直接做 DGF 求逆，以 \Theta(\sqrt n\log n) 的时间花销将其求出。

综上，我们可以在 \Theta(n\log \log n) 时间内求解 1/(2-G)。这个方法同样也适用于求 1/H，其中 H 的乘法可以快速计算。

相关题目：P2025 Dirichlet 半在线卷积。

广义求逆（等比 DGF 求和）

给定 DGF F 和正整数 n,m，我们应该如何快速计算 1+F+F^2+\cdots+F^m 的前 n 项呢？

熟悉等比数列求和的读者会说，1+F+\cdots+F^m=(F^{m+1}-1)/(F-1)，然后采用 LOJ #6713 的方法就行了。不过我们会发现，如果 F[1]=1 的话，F-1 就不可逆，最后会卡在求逆上。而且，因为 DGF 本质是多元多项式，我们不能像一元多项式一样分子分母除个 x 再求逆。那么，我们有没有什么不用倍增的快速做法呢？

有的兄弟，有的。这就是广义求逆。

适当泛化问题，给定 DGF F,G，若已知存在 H 使得 F=GH，且 F[1]=G[1]=0,G[2]\neq 0，那要如何求出 H？

根据 Dirichlet 卷积的式子：

稍微整理一下就是：

依此式递推即可，时间复杂度为 \Theta(n\log n)。

这个做法需要得知 F,G 的前 2n 项，才能算出 H 的前 n 项。看似浪费了很多已知信息，实则我们已经利用了所有信息。因为想要唯一确定 h(n)，就必须知道 f(2n) 是多少。即便是递推式中未出现的 f 的奇数值也很重要，它可以告诉我们这样的 H 是否存在。不过，如果我们已知 H 是存在的，那奇数值确实没啥用（

相关题目：P1998 DGF 等比求和。

后话

感谢 Naszt 和 wkywkywky 对本文内容的贡献。

最初由 Naszt 提出的题目意外地成为契机，让我首次意识到 DGF 牛顿迭代的可能。而在后续的交流中，那些看似零散的讨论碎片，竟一步步拼出了清晰的认知图景。在讨论过程中，Naszt 也提出过另一种截然不同的 \Theta(n\log n) 牛迭算法，可惜囿于算法常数而不够实用。非常值得一提的是，广义求逆的做法就是他提出的。

而 wkywkywky 提出的半在线卷积问题本身已然是块璞玉。最初的探索受数论分治结构的复杂所困，直到无意中瞥见牛顿迭代的晨光穿透云层，这一切才豁然开朗。这场看似偶然的关联，实则昭示着牛顿迭代的精妙普适性。

这段研究的历程长达三个多月，而这所有的内容仅仅始于一个有趣的问题。或许这正是理论研究的迷人之处：某个问题的提出本身，就可能成为照亮新维度的明灯。