Dirichlet series 与算术函数前缀和

Posted on In Disqus: Word count in article: 2k Reading time ≈ 7 mins.

\(\newcommand{\floor}[1]{\left\lfloor #1 \right\rfloor}\)这里我讲讲如何用 Dirichlet series 求一些算术函数的前缀和。这个算法我是在 Project Euler 的论坛里看到的，Lucy_Hedgehog 的帖子。但是他（她？）写的比较粗略，我就稍微写详细一点吧……

Preliminary

这里我们定义几个 notation，并交代几个基础性质。名字什么的我都是乱翻译的，可能和已有的名字不同，如果看到了请告知。

算术函数（Arithmetic function） 若函数 \(f\) 定义域为 \(\mathbb{N}_+\)，值域为 \(\mathbb{C}\)，则称其为算术函数。定义为其前缀和函数 \(S_f(n)_= \sum\limits_{x=1}^n f(x)\)。

我们的问题是给定 \(N\)，如何快速求 \(S_f(N)\)，特别地，我们希望时间复杂度是 sublinear 的。

积性函数（Multiplicative function） 若算术函数 \(f\) 满足：如果 \((n, m) = 1\)，则有 \(f(n) f(m) = f(nm)\)，则称 \(f\) 为积性函数。易知所有积性函数 \(f\) 满足 \(f(1) = 1\)。不过这篇文章没有用到任何积性函数的性质……

Dirichlet series 给定算术函数 \(f\)，定义其 Dirichlet series 为 \[ D_f(s) = \sum_{n=1}^\infty \frac{f(n)}{n^s}. \] 举几个例子： \[ \quad D_\boldsymbol{1}(s) = \sum_{n=1}^\infty \frac{1}{n^s} = \zeta(s), \quad D_\text{Id}(s) = \sum_{n=1}^\infty \frac{n}{n^{-s}} = \zeta(s - 1). \] Dirichlet 卷积（convolution） 定义两个算术函数 \(f, g\) 的 Dirichlet convolution 为 \[ (f * g)(n) = \sum_{d | n} f(d) g(n/d). \] 容易验证以下事实：

  1. Dirichlet 卷积满足交换律、结合律。
  2. 令 \(\varepsilon(n) = [n = 1]\)，则对于任意算术函数 \(f\)，均有 \(\varepsilon * f = f\)。

较为经典的几个例子有 \[ (\phi * \boldsymbol{1})(n) = \sum_{d | n} \phi(d) = n = \text{Id}(n). \] 以及 \[ (\mu * \boldsymbol{1})(n) = \sum_{d|n} \mu(d) = [n = 1] = \varepsilon(n). \] Dirichlet 逆元（inverse） 若 \(f * g = \varepsilon\)，则称 \(g\) 为 \(f\) 的 Dirichlet 逆元，记做 \(g = f^{-1}\)。这里提供了很多逆元的例子。上面提到 \(\mu\) 的 Dirichlet 逆元就是 \(\boldsymbol{1}\)。Möbius 反演可以简单理解为 \[ f * \boldsymbol{1} = g, \quad \Leftrightarrow \quad f = g * \mu. \] 注意 Dirichlet 逆元不一定存在。

Dirichlet convolution 与 Dirichlet series 令 \(f, g\) 为两个算术函数，有 \[ D_{f * g}(s) = \sum_{n=1}^\infty \frac{1}{n^s} \sum_{d | n} f(d) g(n/d) = \sum_{a, b \geq 1} \frac{f(a) g(b)}{(ab)^s} = \sum_{a=1}^{\infty} \frac{f(a)}{a^s} \sum_{b=1}^\infty \frac{g(b)}{b^s} = D_f(s) D_g(s). \] 即 \(f * g\) 的 Dirichlet series 为 \(f\) 的 Dirichlet series 乘上 \(g\) 的 Dirichlet series。

算法 high-level idea

我们先固定一个 \(N \in \mathbb{N}\)。对于一个算术函数 \(f\)，我们 somehow 维护一个数据结构来表达 \(D_f(s)\)，并且这个数据结构支持以下操作：

  1. 给定 \(D_f(s)\) ，能快速计算 \(S_f(N)\)。注意这个数据结构不需要对于任何一个 \(n\) 都给出 \(S_f(n)\)。
  2. 对于一些较为“简单”的算术函数 \(f\)，能够快速初始化这个数据结构。
  3. 给定 \(D_f(s)\) 和 \(D_g(s)\)，能快速计算 \(D_{f + g}(s), D_{f - g}(s), D_{f * g}(s), D_{f / g}(s)\)。

那么我们就可以把目标算术函数 \(f\) 的 Dirichlet series 用若干个较为“简单”的算术函数的 Dirichlet series 的四则运算表达出来，然后直接查询 \(S_f(N)\) 即可。

算法详情

这个数据结构非常简单：令 \(w\) 为一待取常数，我们记录所有的 \(\{f(n)\}_{n=1}^w\) 以及所有的 \(\left\{S_f(\floor{\frac{N}{i}}) \right\}_{i=1}^{\floor{N/w}}\)。这样查询非常简单，因为我们已经记录了 \(S_f(N)\)。接下来我们介绍如何支持第三个操作。由于 \(D_{f+g}(s), D_{f-g}(s)\) 过于 trivial，我就懒得写了……

计算 \(D_{f*g}(s)\)

不妨令 \(h = f * g\)。对于 \(\{h(n)\}_{n=1}^w\)，直接用定义暴力就好了，复杂度是 \(O(w \log w)\)。

对于 \(\left\{S_h(\floor{\frac{N}{i}}) \right\}_{i=1}^{\floor{N/w}}\) 这一块，其实也挺好整的：对于任意 \(x \in \mathbb{N}_+\)， \[ \begin{equation} S_h \left(x \right) = \sum_{n \leq x} h(n) = \sum_{a, b: ab \leq x} f(a) g(b) = \sum_{a \leq x} f(a) \sum_{b \leq x/a} g(b) = \sum_{a \leq x} f(a) S_g\left( \floor{\frac{x}{a}} \right). \label{eq:S_c} \end{equation} \] 由于 \(x = \floor{N/i}\)，所以实际上我们已经有 \(S_g(\floor{x/a})\) 了，且所有不同的 \(\floor{x/a}\) 只有 \(O(\sqrt{x})\) 个。故我们使用经典 trick，可以枚举出所有不同的 \(\floor{x/a}\)，对一个区间内的 \(x\) 利用 \(S_f\) 一次计算完毕。这样计算一个 \(S_h(\floor{N/i})\) 的复杂度为 \(O(\sqrt{N/i})\)，计算出所有 \(S_h(\floor{N/i})\) 的复杂度为： \[ \sum_{i=1}^{\floor{N/w}} \sqrt{N/i} \approx \int_{i=1}^{N/w} \sqrt{N/i} \mathrm d i = O(N / \sqrt{w}). \] 这里直接用 Dirichlet hyperbola method 也行，复杂度不变： \[ S_h(x) = \sum_{a \leq \sqrt{x}} f(a) S_g(\floor{x / a}) + \sum_{a \leq \sqrt{x}} g(a) S_f(\floor{x/a}) - S_f(\floor{\sqrt{x}}) S_g(\floor{\sqrt{x}}). \] 整体复杂度为 \(O(w \log w + N / \sqrt{w})\)。

计算 \(D_{h/f}(s)\)

过程和 \(D_{f*g}(s)\) 基本上类似。

令 \(g = h / f\)。对于 \(\{g(n)\}_{n=1}^w\)，可得 \[ \sum_{d|n} f(d) g(n/d) = h(n) \quad \Rightarrow \quad f(1) g(n) = h(n) - \sum_{d | n, d > 1} f(d) g(n/d). \] 故依次计算 \(\{g(n)\}_{n=1}^w\) 即可。复杂度依然是 \(O(w \log w)\)。写法上可能需要注意。

对于第二部分，由 \(\eqref{eq:S_c}\) 可以直接推出 \[ f(1) S_g \left(x \right) = S_h \left(x \right) - \sum_{1 < a \leq x} f(a) S_g\left( \floor{\frac{x}{a}} \right), \] 复杂度分析类似乘法，也是 \(O(w \log w + N / \sqrt{w})\)。

从这个算法也可以看出，\(f^{-1}\) 存在当且仅当 \(f(1) \neq 0\)。

取 \(w\)

很显然，我们只要让 \(w \log w = N / \sqrt{w}\)，即 \(w = \Theta((N / \log N)^{2/3})\)，就可以使得整体复杂度取得最小，为 \(O(N^{2/3} \log^{1/3} N)\)。

这个算法里面我没有仔细想 \(O(w \log w)\) 这部分是否可以优化到 \(O(w)\)，如果可以的话复杂度可以再降一个 \(O(\log ^{1/3} N)\)。比起这篇文章来说，复杂度降低了 \(O(N^{1/12}/\log^{1/3} N)\)，可喜可贺，可喜可贺。不过我感觉那篇文章的参数选的有问题，他用的 \(w = N^{1/2}\)……

“简单”的算术函数

对于一些函数，我们能够直接求出它的 Dirichlet series 表示形式，这些函数就可以作为算法的起点。常用的有以下几个例子：

  1. \(\text{Id}_k(n) = n^k\) 的 Dirichlet series 为 \(\zeta(s - k)\)，可以直接初始化该数据结构。\(k\) 大的时候可能需要 Faulhaber 公式。
  2. 若 \(f\) 只在少数几个输入上非 0，则也可以直接初始化该数据结构。

简单例题

\[ D_\phi(s) = \frac{D_{\text{Id}}(s)}{D_\boldsymbol{1}(s)} = \frac{\zeta(s - 1)}{\zeta(s)}. \] 于是 \(\phi\) 的前缀和可以快速搞定。

\[ \mu * \boldsymbol{1} = \varepsilon, \quad \Rightarrow \quad D_\mu(s) = \frac{1}{\zeta(s)}， \] 于是 Mertens function 也可以在 \(O(N^{2/3} \log^{1/3} N)\) 内搞定。注意到我们在算 \(D_\mu(s)\) 的时候，\(\{\mu(n)\}_{n=1}^w\) 这一块可以直接线性筛出来而不是暴力用除法搞，这样这部分的复杂度可以从 \(O(w \log w)\) 降到 \(O(w)\)，取 \(w = N^{2/3}\) 就可以把整体时间复杂度降到 \(O(N^{2/3})\)，和经典方法一样了。

这类算法最大的问题是如何用简单的算术函数通过给定的运算表达出目标算术函数，我试了几个题，表示有点难……主要是 Dirichlet inverse 根本看不出有啥性质……

对 Project Euler 有兴趣的同学可以尝试 PE 351/415/428/447/448/608/715，还可以看一下 Lucy_Hedgehog 和 MuthuVeerappanR 在 428 下的 post，以及 715 下 Jaehyun 的 post。我猜 MuthuVeerappanR 就是刚才提到的那篇 blog 的作者。

代码的话我找时间 push 到 github 上去吧……

这文章写的比之前的《解析方法数质数：从入门到入土》轻松多了……