DGF 的计算理论：经典方法和牛顿迭代

· · 算法·理论

DGF 的计算理论：经典方法和牛顿迭代

观前提醒

这篇文章一定程度上是 飞雨烟雁 的 这篇文章 的详解，
你可以先阅读完那篇文章再看本文以做出更深刻的理解。

我们的研究成果是对我出的一题的探究，
这题本来要作为 我的比赛 的 E 题，不过由于难度过大，就干脆将其公开共大家一起探讨。

简介

在普通多项式科技中，我们对于一个能写成关于 OGF/EGF 的方程，可以使用牛顿迭代解决。

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

在 DGF 中，有一类筛法问题，试求积性函数 F 的某项的前缀和，
可以用杜教筛等方法：预处理前 n^{2/3} 项，再进行杜教筛。

本文我们解决的是对于 DGF 的方程，试求多项式 F 的某项的前缀和。

我用比较自然的推导给出了一种较为经典的方法，
而飞雨烟雁再次发现了 DGF 的牛顿迭代，这种方法的编码难度也低很多。

如果用杜教筛作为基本的卷积方法，那么对于一个一般的低次方程：
递推计算部分的时间复杂度为：\Theta(T\log T)
块筛计算部分的时间复杂度为：\Theta(T^{-1/2}n)
即为 \Theta(n^{2/3}\log n) 的总时间复杂度。

如果用 zak 筛作为基本的卷积方法，时间复杂度理论上可以优化到 \tilde O(\sqrt n) 的时间复杂度。

题目

本文依然用我们研究的题目作为具体的例子进行讲解，
但在此我就不具体的给出题目了，这是该题的形式化描述：

给定 n,m,f_1=1，f_n 的递推式为：

f_n=\sigma_1(n) +\sum_{\substack{d_1=n\\[3pt]d_1<n}}f_{d_1} +\sum_{\substack{d_1d_2=n\\[3pt]d_1,d_2<n}}f_{d_1}f_{d_2} +\dots +\sum_{\substack{d_1d_2\dots d_m=n\\[3pt]d_1,d_2,\dots,d_m<n}}f_{d_1}f_{d_2}\dots f_{d_m}

求：

ans=\sum_{i=1}^nf_i

定义 f_n 的 DGF：

F(x)=\sum_{n\ge1}\frac{f_n}{n^x}

可以得到：

\left(\frac{m^2}{2}+\frac{m}{2}+1\right)F(x)=\frac{m^2}{2}-\frac{m}{2}+\zeta(x)\zeta(x-1)+F(x)+F^2(x)+\dots+F^m(x)

我个人认为这个方程是比较有代表性的，难以找到什么逃课方法。

经典方法

\begin{aligned} \left(\frac{m^2}{2}+\frac{m}{2}+1\right)F(x)&=\frac{m^2}{2}-\frac{m}{2}+\zeta(x)\zeta(x-1)+F(x)+F^2(x)+\dots+F^m(x)\\ \left(\frac{m^2}{2}+\frac{m}{2}+1\right)F(x)-\zeta(x)\zeta(x-1)-\frac{m^2}{2}+\frac{m}{2}+1&=\frac{F^{m+1}(x)-1}{F(x)-1}\\ F^{m+1}(x)-\left(\frac{m^2}{2}+\frac{m}{2}+1\right)F^2(x)+m^2F(x)&=\zeta(x)\zeta(x-1)+\frac{m^2}{2}-\frac{m}{2}-\zeta(x)\zeta(x-1)F(x) \end{aligned}

递推计算

定义：

f_n^m=[n^{-x}]F^m(x)=\sum_{d_1d_2\dots d_m=n}f_{d_1}f_{d_2}\dots f_{d_m}

生成函数方程可化为：

f_n^{m+1}-\left(\frac{m^2}{2}+\frac{m}{2}+1\right)f_n^2+(m^2+1)f_n=\sigma_1(n)-\sum_{\substack{d\mid n\\[3pt]d<n}}\sigma_1(d)f_d

根据卷积原理:

f_n^{a+b}-f_n^a-f_n^b=\sum_{\substack{d\mid n\\[3pt]1<d<n}}f_d^af_{\frac{n}{d}}^b

我们以 n 递增的顺序来计算，这样对于任意一个时刻，右式都为常数。
一种方法是引入 m 个由卷积原理得到的式子，得到一个 m+1 元方程组，再进行消元。
但我们发现系数矩阵是不和 n 有关的，我们可以提前消好元，每次直接计算即可。

但我在前面之所以乘上 (F(x)-1) 的原因之一就在此：
由于任意两相邻数字可被写成别的两相邻数字之和。
按照这样倍增的思想，我们仅需引入 \Theta(2\log m) 个方程和未知数即可完成消元。

块筛计算

定义：

s_n^m=\sum_{i=1}^nf_i^m

生成函数方程可化为：

s_n^{m+1}-\left(\frac{m^2}{2}+\frac{m}{2}+1\right)s_n^2+(m^2+1)s_n=\sum_{d=1}^n\sigma_1(d)+\frac{m^2}{2}-\frac{m}{2}-\sum_{i=2}^n\sigma_1(i)s_{\lfloor n/i\rfloor}

根据杜教筛原理：

s_n^{a+b}-s_n^a-s_n^b=\sum_{i=2}^{\lfloor n/2\rfloor}\Delta s_i^as_{\lfloor n/i\rfloor}^b-s_{\lfloor n/2\rfloor}^a

按照一样的方法消元即可。

BUG-线性相关

由于这题的生成函数的特殊性， 这个方法对于这题还是有个小 BUG。

注意到我们无法解方程，方程是线性相关的。
本质上是因为方程无法提供关于第 n 项的信息：

\begin{aligned} f_n^{m+1}-\left(\frac{m^2}{2}+\frac{m}{2}+1\right)f_n^2+(m^2+1)f_n&=\sigma_1(n)-\sum_{\substack{d\mid n\\[3pt]d<n}}f_d\\ \sum_{\substack{d_1d_2\dots d_{m+1}=n\\[3pt]d_1,d_2,\dots,d_{m+1}<n}}f_{d_1}f_{d_2}\dots f_{d_{m+1}} -\left(\frac{m^2}{2}+\frac{m}{2}+1\right)\sum_{\substack{d_1d_2=n\\[3pt]d_1,d_2<n}}f_{d_1}f_{d_2} &=\sigma_1(n)-\sum_{\substack{d\mid n\\[3pt]d<n}}f_d \end{aligned}

这个方程中其实没有一项和 f_n 有关。

但是通过消元我们可以得到一个不含 f_n^x 的恒等式：

\sum_{i=0}^v\left(k_i\sum_{\substack{d\mid n\\[3pt]1<d<n}}f_d^{l_i}f_{\frac{n}{d}}^{r_i}\right)=\sigma_1(n)-\sum_{\substack{d\mid n\\[3pt]d<n}}f_d

于是我们可以通过 2n 凑出 n 这一项：

\begin{aligned} \sum_{i=0}^v\left(k_i\sum_{\substack{d\mid 2n\\[3pt]1<d<2n}}f_d^{l_i}f_{\frac{2n}{d}}^{r_i}\right)&=\sigma_1(n)-\sum_{\substack{d\mid 2n\\[3pt]d<2n}}f_d\\ \sum_{i=0}^v\left(k_i\sum_{\substack{d\mid 2n\\[3pt]2<d<n}}f_d^{l_i}f_{\frac{2n}{d}}^{r_i} +k_i(f_2^{l_i}f_n^{r_i}+f_n^{l_i}f_2^{r_i}) \right)&=\sigma_1(n)-\sum_{\substack{d\mid 2n\\[3pt]d<n}}f_d-f_n\\ \sum_{i=0}^v\left(k_if_n^{r_i}+k_if_n^{l_i}\right)+f_n&=\sigma_1(n)-\sum_{\substack{d\mid 2n\\[3pt]d<n}}f_d-\sum_{i=0}^v\left(k_i\sum_{\substack{d\mid 2n\\[3pt]2<d<n}}f_d^{l_i}f_{\frac{2n}{d}}^{r_i}\right) \end{aligned}

对于这个方程就可以正常消元了。

杜教筛部分也同理：

\begin{aligned} \sum_{i=0}^v\left(k_i\sum_{i=2}^{\lfloor 2n/2\rfloor}\Delta s_i^{l_i}s_{\lfloor 2n/i\rfloor}^{r_i}-k_is_{\lfloor 2n/2\rfloor}^{l_i}\right)&=\sum_{d=1}^{2n}\sigma_1(d)+\frac{m^2}{2}-\frac{m}{2}-\sum_{i=2}^ns_{\lfloor 2n/i\rfloor}\\ \sum_{i=0}^v\left(k_i\sum_{i=3}^{\lfloor 2n/3\rfloor}\Delta s_i^{l_i}s_{\lfloor 2n/i\rfloor}^{r_i}+k_i(s_n^{r_i}+s_{\lfloor n\rfloor}^{l_i}-2s_{\lfloor 2n/3\rfloor}^{l_i})\right)&=\sum_{d=1}^{2n}\sigma_1(d)+\frac{m^2}{2}-\frac{m}{2}-\sum_{i=3}^ns_{\lfloor 2n/i\rfloor}-s_n\\ \sum_{i=0}^v\left(k_is_n^{r_i}+k_is_n^{l_i}\right)+s_n&=\sum_{d=1}^{2n}\sigma_1(d)+\frac{m^2}{2}-\frac{m}{2}-\sum_{i=3}^ns_{\lfloor 2n/i\rfloor}-\sum_{i=0}^v\left(k_i\sum_{i=3}^{\lfloor 2n/3\rfloor}\Delta s_i^{l_i}s_{\lfloor 2n/i\rfloor}^{r_i}+2k_is_{\lfloor 2n/3\rfloor}^{l_i}\right)\\ \end{aligned}

时间复杂度

有 O(\log m) 个方程，相当于做 O(\log m) 次杜教筛。
时间复杂度是 O(B\log n\log m+B^{-1/2}n\log m)，
取 B=n^{2/3} 时较优，即为 O(n^{2/3}\log n\log m)。

牛顿迭代

递推计算

DGF 牛顿迭代

仿照贝尔级数的思想，我们可以把 p^{-x} 当做一个单项式 x_{p}。
这便是多元多项式问题了。

G(x):=-\left(\frac{m^2}{2}+\frac{m}{2}+1\right)x+\frac{m^2}{2}-\frac{m}{2}+\sigma_1+x+x^2+\dots+x^m

那么我们可以根据多项式牛顿迭代求解出 F(x)：

\mathbf x:=(x_1,x_2,\cdots,x_m) G(F(\mathbf x))\equiv 0\pmod {\mathbf x^n}

但是实际计算时我们根本不需要用多元多项式，
我们只是借助多元多项式这个框架来说明牛顿迭代的适用性。

那么由于指数最高时也不超过 \log_2n，而牛顿迭代的精度是翻倍的，所以我们仅用 \log\log n 次牛顿迭代即可。

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

但实际上由于我们做了非常多没意义的运算，优化一下转移顺序，不做这些没意义的运算就可以做到 \Theta(n\log n) 的时间复杂度了。
这个方法和经典做法也是及其类似的，可惜囿于算法常数而不够实用。

但是有一种工程上更方便的方法：
我们可以按照一元多项式迭代的方法进行计算，飞雨烟雁将其称为非齐次牛顿迭代。
这样的时间复杂度也优化到了同样的 T(n)=T(\sqrt n)+\Theta(n\log n)=\Theta(n\log n)。

函数复合

但现在就真的可以做到 \Theta(n\log n) 的时间复杂度了吗？

牛顿迭代的具体步骤是：

F(\mathbf x)\equiv F_0(\mathbf x)-\frac{G(F_0(\mathbf x))}{G'(F_0(\mathbf x))}\pmod {\mathbf x^{2n}} \begin{aligned} G(x)&=-\left(\frac{m^2}{2}+\frac{m}{2}+1\right)x+\frac{m^2}{2}-\frac{m}{2}+\sigma_1+x+x^2+\dots+x^m\\ G'(x)&=-\left(\frac{m^2}{2}+\frac{m}{2}+1\right)+1+2x+3x^2+\dots+mx^{m-1}\\ \end{aligned}

可以发现我们需要计算 x+x^2+\dots+x^m 和 1+2x+3x^2+\dots+mx^{m-1}，
而 DGF 卷积是 \Theta(n\log n) 的，所以这里 G\circ F 要 \Theta(nm\log n) 的时间复杂度。

有一种方法是：

\begin{aligned} G(x)&=-\left(\frac{m^2}{2}+\frac{m}{2}+1\right)x+\frac{m^2}{2}-\frac{m}{2}+\sigma_1-1+\frac{x^{m+1}-1}{x-1}\\ G'(x)&=-\left(\frac{m^2}{2}+\frac{m}{2}+1\right)+\frac{mx^{m+1}-(m+1)x^m+1}{(x-1)^2}\\ \end{aligned}

分子可用 DGF 快速幂 \Theta(n\log n) 计算。

然而 [1^{-x}]G=1，即 [1^{-x}](G-1)=0，我们无法对其求逆。
但可以回到本质上来看这个除法，我们将其称为广义求逆。
可以发现广义求逆的方法和旧做法解决线性相关的方法本质相同。
时间复杂度是 \Theta(n\log n) 的。

不过后面会发现这里并不是时间复杂度的大头所在，所以没必要用这么麻烦的方法。

另一种工程上更简单方法是用倍增，这就和经典方法比较类似了：

\begin{aligned} F_p&=F^m\\ F_s&=F+F^2+\dots+F^m\\ F_t&=1+2F+3F^2+\dots+mF^{m-1}\\ \end{aligned}

类似于快速幂，对于 2m 可以这么计算：

\begin{aligned} F_p&\leftarrow F_p^2\\ F_s&\leftarrow F_s(F_p+1)\\ F_t&\leftarrow F_t(F_p+1)+mF_p(F_s-F_p+1)\\ \end{aligned}

时间复杂度是 \Theta(n\log m\log n) 的。

块筛计算

我们先预处理出前 \sqrt n 项，
根据非齐次牛顿迭代的原理，只用再迭代一次即可：

F(\mathbf x)\equiv F_{\le\sqrt n}(\mathbf x)-\frac{G(F_{\le\sqrt n}(\mathbf x))}{G'(F_{\le\sqrt n}(\mathbf x))}\pmod {\mathbf x^{n}}

我们想求 F 的前缀和 s_n，求出等式右边的块筛处的值即可，
※ 这里和下文中块筛都指的是 n 处的块筛。
若 H=F*G 那么可以根据 F 和 G 的块筛处的值求出 H 块筛处的值。

同样也可以根据 某函数块筛处的值 求出 某函数的逆的块筛处的值。

这样我们预处理出前 n^{2/3} 项，用朴素的杜教筛可以做到 O(n^{2/3}\log n\log m)。

zak 筛简介

根据 周康阳 -《关于积性函数求和问题的一些进展》所述，
可以在 \tilde O(\sqrt n) 的时间复杂度完成这类块筛卷积，块筛求逆。

这里计算 G\circ F 时无法使用第一种方法，但第二种倍增的方法还是能用的。

朴素 zak 筛

我这里主要是复读 周康阳的论文 和 相关 blog：

我们要解决的是 h(z)=\sum\limits_{xy=z}f(x)g(y) 在所有 \left\lfloor\frac nt\right\rfloor 上的前缀和。

用朴素的杜教筛可以做到 O(n^{2/3})，
而 zak 筛的主要思想我觉得是考虑主动转移 + FFT 优化：

当 x>\sqrt n 或 y>\sqrt n 时

考虑枚举 y>\sqrt n。
※ x>\sqrt n 时是对称的。

S_{h}\left(\left\lfloor\frac nt\right\rfloor\right)=\sum\limits_{i=1}f(i)S_g\left(\left\lfloor\frac n{it}\right\rfloor\right)

考虑枚举 i,t，因为 it\le\sqrt n，所以时间复杂度是 O(\sqrt n\log n) 的。

当 x<\sqrt n 且 y<\sqrt n 时

不妨同时开对数以使用多项式优化：

h(z)=\sum_{\ln(xy)=\ln(z)}f(x)g(y) h(z)=\left[x^{\ln(z)}\right]\left(\sum_{i=1}^{\sqrt n}f(i)x^{\ln(i)}\right)\left(\sum_{i=1}^{\sqrt n}g(i)x^{\ln(i)}\right)

不过 FFT 的的指数要是整数，我们可以取个估计：
我们认为 ln(x)+\ln(y)\le\ln\left(\frac nt\right) 当且仅当

这样值域就是 $\lceil B\ln(n)\rceil$ 的，用 FFT 可以做到 $\Theta(B\log^2n)$。 既然是估计，那么就有误差，产生误差的充要要条件是： $$\lceil Bln(x)\rceil+\lceil B\ln(y)\rceil>\left\lceil B\ln\left(\frac nt\right)\right\rceil$$ 可以算出一个必要条件是： $$xyt\in\left(n-\frac{2n}B,n\right]$$ 所以我们可以使用区间筛 $\Theta(\frac nB\log n)$ 质因子分解后 dfs 搜索确实产生误差的部分进行更新。 又因为 $\sum_{ijk\le n}1=\Theta(n\log^2n)$，所以暴力 dfs 枚举的时间复杂度是 $\Theta(\frac nB\log^2n)$ 的。 我们取 $B=\sqrt n$ 有 $\Theta(\sqrt n\log^2n)$ 的时间复杂度。 ### 块筛求逆 我们需要计算 $H=\frac FG$， 有 $f(n)=\sum\limits_{ij=n}h(i)g(j)$。 记 $H_{\le\sqrt n}$ 是 $H$ 只保留 $\le\sqrt n$ 位置上的值的函数。 考虑枚举 $j\le\sqrt n$，可以 $\Theta(\sqrt n\log n)$ 预处理 $H_{\le\sqrt n}$。 那么有 $F=(H_{\le\sqrt n}+H_{>\sqrt n})G$，可以通过一次块筛卷积得到 $H_{\le\sqrt n}G$。 记 $A=F-H_{\le\sqrt n}G$，有： $$S_{a}\left(\left\lfloor\frac nt\right\rfloor\right)=\sum\limits_{i=1}g(i)S_{h_{>\sqrt n}}\left(\left\lfloor\frac n{it}\right\rfloor\right)$$ 考虑枚举 $i,t$，时间复杂度是 $O(\sqrt n\log n)$ 的。 ## 总结 我们把这些方法综合起来就可以有 DGF 多项式 / 块筛 全家桶了: 普通多项式全家桶： > NTT > operator +-*/ > der,inte,inv,exp,ln > Newton_Iteration DGF 多项式全家桶： > operator +-*/ 「这里的除法可以写广义求逆」 > der,inte 「我们并不在意真实值，形式上有相同的性质即可」 > exp,ln > Newton_Iteration 块筛全家桶： > operator +- > operator * 「zak 筛理论更快，不过实际表现更差，可以写杜教筛」 > operator / 「可以查看上文的块筛求逆」 这样对于一般的，可以用杜教筛解决的题目都可以套模板了（但是常数巨大）， 如果块筛全家桶中再实现 zak 论文中的 $ln,exp$ 用途就更广了。 ## 后话 感谢 [飞雨烟雁](https://www.luogu.com.cn/user/375984) 对本文内容的贡献。 我个人觉得 DGF 技术的进一步飞跃是比较困难的， 毕竟狄利克雷卷积并不常见，题目也需要精妙构造。 在 OGF 中我们有简单的移位算子 —— 乘 $x$。 而 DGF 一旦出现加法相关的东西就及其复杂了。