解析方法数质数：从入门到入土（下）

Posted on In Disqus: Word count in article: 6k Reading time ≈ 22 mins.

\(\newcommand{\eps}{\epsilon} \newcommand{\tO}{\widetilde O} \newcommand{\d}{\mathrm d} \newcommand{\hphi}{\widehat \phi} \newcommand{\hPhi}{\widehat \Phi} \newcommand{\sinc}{\text{sinc}}\)上次我们讲了 [LO]，[Galway] 和 [Platt] 的大致思路，似乎一切都已经明了，只剩下最后一个难点了：\(\zeta\) 函数怎么求？这里我们就来介绍一下我 survey 到的几种求 \(\zeta(s)\) 的方法。

这里我们先约定：文中的所有 \(\sigma\) 可以被替换为 \(\Re(s)\)，\(t\) 可以被替换为 \(\Im(s)\)。

我们再对 \(\zeta(s)\) 的需求细分一下：在 [LO] 和 [Galway] 中，我们真的需要 \(\zeta(s)\) 的值，因为最后的积分需要，而且这里的 \(s\) 还不满足 \(\sigma = \frac{1}{2}\)。而在 [Platt] 中，我们并不关心 \(\zeta(s)\) 的值，我们只关心 \(\zeta(s)\) 的非平凡零点在哪里。

想必大家都听说过 Riemann hypothesis 的名字。它是数学界最著名的 hypothesis 之一，说的就是：\(\zeta(s)\) 的所有非平凡零点均有 \(\sigma = \frac{1}{2}\)。于是大家便称 \(\sigma = \frac{1}{2}\) 这条线为 critical line。我们在接下来的分析中也会经常看到 critical line 上的 \(s\) 会有特殊的性质。

Euler-Maclaurin 公式

TL;DR：用 \(\eqref{eq:euler-maclaurin-error-term}\) 来计算一个合适的 \(m\)，然后用 \(\eqref{eq:zeta-euler-maclaurin}\) 来计算 \(\zeta(s)\)。

当然最古老是还是 Euler-Maclaurin 公式。在二十世纪三十年代之前一直都是用这个公式来计算 \(\zeta(s)\) 的。我们再来复习一下 Euler-Maclaurin 公式：

[Euler-Maclaurin Formula] 令 \(a< b\) 为两自然数，对于任意正整数 \(m\)，\(f\) 在 \([a, b]\) 上连续，则有 \[ \int_a^b f(x) \d x = \sum_{i=a}^b f(x) - \frac{1}{2}(f(a) + f(b)) - \sum_{k=1}^{m-1} \frac{B_{2k}}{(2k)!}(f^{(2k-1)}(b) - f^{(2k-1)}(a)) + R_{2m}, \] 其中 \(B_{2k}\) 为 Bernoulli 数，\(R_{2m}\) 为余项。

我们不妨把取 \(f(x) = x^{-s}\)，有 \(f^{(2k-1)}(x) = -x^{1-2k-s} \prod\limits_{j=0}^{2k-2} (s + j)\)，带入即有 \[ \int_a^b x^{-s} \d x = \sum_{i=a}^b x^{-s} - \frac{a^{-s} + b^{-s}}{2} + \sum_{k=1}^{m-1} \frac{B_{2k}}{(2k)!} \left( b^{1-2k - s} - a^{1-2k - s} \right) \prod_{j=0}^{2k-2} (s+j) + R_{2m}, \] 然后我们取 \(b\) 的极限到 \(\infty\)，有 \[ \begin{equation} \zeta(s) = \sum_{i=1}^\infty x^{-s} = \sum_{i=1}^{a-1} x^{-s} + \int_a^\infty x^{-s} \d s + \frac{a^{-s}}{2} + \sum_{k=1}^{m-1} \frac{B_{2k}}{(2k)!} a^{1-2k-s} \prod_{j=0}^{2k-2} (s+j) + R_{2m}. \label{eq:zeta-euler-maclaurin} \end{equation} \] 中间这个积分非常显然，我就懒得化开了。不妨令 \(T_{a, k}(s) := \frac{B_{2k}}{(2k)!} a^{1-2k-s} \prod\limits_{j=0}^{2k-2} (s+j)\)，在 \(\sigma > -2m + 1\) 的情况下，[Eq 1.3, OS] 给了 \(R_{2m}\) 的一个 bound： \[ \begin{equation} |R_{2m}| \leq \left| \frac{s + 2m - 1}{\sigma + 2m - 1} T_{a, m}(s) \right|. \label{eq:euler-maclaurin-error-term} \end{equation} \] 故我们需要选择合理的 \(m\) 和 \(a\) 使得 \(R_{2m}\) 小。当 \(\Re(s) = \frac{1}{2}\) 时，可以选择 \(a \approx t / (2\pi)\)，\(m \approx \sqrt{n}\)，不过我也没算这样是否合理……但是从这里可以看出：Euler-Maclaurin 公式的复杂度是 \(\tO(t)\) 的，在 \(t\) 大的情况非常菜，所以我们还需要改进……

Riemann-Siegel 公式

注意到我之前说二十世纪三十年代之前都是用 Euler-Maclaurin 公式来计算 \(\zeta(s)\) 的，那么是什么打破了这一局面呢？没错就是 Siegel 在整理 Riemann 在 1850 年代的手稿时发现了这个公式，于是被称为 Riemann-Siegel 公式。

首先我们来介绍一个 Riemann-Siegel 公式。它有很多变种，我们先说其中一个：[Eq 7.2, Galway], [Eq 1.1, de Reyna] 对于任意 \(s\)，有 \[ \zeta(s) = \mathcal{R}(s) + \chi(s) \overline{\mathcal{R}(1 - \bar s)}, \quad \mathcal{R}(s) := \int\limits_{0 \swarrow 1} \frac{z^{-s} e^{i \pi z^2}}{e^{i\pi z} - e^{-i\pi z}} \d z, \] 其中 \(0 \swarrow 1\) 描述 \(\frac{1}{2} - e^{i \pi /4} \alpha\) 这条路径，\(-\infty < \alpha < \infty\), \(\chi(s)\) 为 \[ \chi(s) := \pi^{s - 1/2} \frac{\Gamma((1 - s) / 2)}{\Gamma(s / 2)}. \] 注意到 \(f(z; s) := \frac{z^{-s} e^{i \pi z^2}}{e^{i \pi z} - e^{-i \pi z}}\) 是一个亚纯函数，极点为所有整数，且在正整数 \(n\) 上的留数为： \[ \text{Res}(f; n) = \lim_{z \to n} f(z; s) (z - n) = z^{-s} e^{i \pi z^2} \lim_{z \to n} \frac{z - n}{e^{i \pi z} - e^{-i \pi z}} = \frac{z^{-s}}{2 \pi i}. \] 故我们可以把 \(0 \swarrow 1\) 这条线移到 \(N \swarrow N+1\) ，由留数定理可得： \[ \mathcal{R}(s) = \int\limits_{0 \swarrow 1} f(z; s) \d z = 2 \pi i \sum_{n=1}^N \text{Res}(f; n) + \int\limits_{N \swarrow N+1} f(z; s) \d z = \sum_{n=1}^N n^{-s} + \int\limits_{N \swarrow N+1} f(z; s) \d z. \] 对于任意非负整数 \(N, M\)，均有 \[ \begin{equation} \zeta(s) = \sum_{n=1}^N n^{-s} + \chi(s) \sum_{n=1}^M n^{s-1} + R_{N, M}(s). \label{eq:RS-1} \end{equation} \] 其中 \(R_{N, M}\) 为余项。我们可以通过分析余项来分析 \(\eqref{eq:RS-1}\) 的准确度。

[Gabcke] 在 Critical Line 上的 \(\zeta(s)\) 算法

TL;DR：用 \(\eqref{eq:RS-zeta}\) 来计算 \(\theta(t)\)，用 \(\eqref{eq:gabcke}\) 来计算 \(Z(t)\)，用 \(\eqref{eq:gabcke-error-term}\) 计算在何处截断 \(\eqref{eq:gabcke}\) 里的级数，最后用 \(\eqref{eq:RS-zeta-Z-theta}\) 来计算 \(\zeta(1/2+it)\)。

我们先写出一个关于 \(\zeta(s)\) 的美妙的对称性质： \[ \begin{equation} \pi^{-s/2} \Gamma(s/2) \zeta(s) = \pi^{-(1-s) / 2} \Gamma((1-s)/2) \zeta(1 - s). \label{eq:functional} \end{equation} \] 这个对称性揭示了 \(\xi(s) := \pi^{-s/2} \Gamma(s/2) \zeta(s)\) 关于 critical line 是对称的，这也就意味着：对于 critical line 上的 \(s\)，有 \[ \overline{\xi(s)} = \xi(\bar s) = \xi(1 - s) = \xi(s), \] 即 \(\xi(s)\) 一定是实数。所以 critical line 上的 \(s\) 有 \(\arg (\pi^{-s/2} \Gamma(s/2)) = -\arg \zeta(s)\)，从而我们可以定义 Riemann-Siegel \(\theta\) 函数： \[ \theta(t) := \arg(\pi^{-(1/2+it)/2} \Gamma((1/2+it)/2)) = \arg \left(\Gamma\left( \frac{1/2 + it}{2} \right) \right) - \frac{1}{2} t \log \pi, \] 来描述 \(\zeta(s)\) 的幅角（的相反数），并以此定义 Riemann-Siegel Z 函数 \[ \begin{equation} Z(t) := e^{i \theta(t)} \zeta\left( 1/2+it \right). \label{eq:RS-zeta-Z-theta} \end{equation} \] 根据前文所述，我们可以知道 \(Z(t)\) 一定是实数。\(Z(t)\) 在研究 \(\zeta(s)\) 在critical line 上的零点是特别有用，因为 \(|Z(t)| = |\zeta(1/2+it)|\)。\(\theta(t)\) 和前文定义的 \(\chi(s)\) 有紧密联系：对 \(\eqref{eq:functional}\) 稍做移项，对于 critical line 上的 \(s\) 有 \[ \chi(s) = \frac{\zeta(s)}{\zeta(1 - s)} = \frac{|Z(t)| \exp(-i \theta(t))}{|Z(-t)| \exp(i\theta(t))} = \exp(-2i \theta(t)). \] [Gabcke] 给出了 \(Z(t)\) 的 Riemann-Siegel 公式：（感受一下即可，背不下来的……劝退开始）

[Theorem 2.1.1, Theorem 2.1.6, Gabcke] \[ \begin{equation} Z(t) = 2 \sum_{n=1}^N \cos(\theta(t) - t \log n) n^{-1/2} + \frac{(-1)^{(N-1)}}{\sqrt{a}} \sum_{n=0}^\infty \frac{C_n(z)}{a^n}, \label{eq:gabcke} \end{equation} \] 其中 \(a := \sqrt{\frac{t}{2\pi}}\)，\(N := \lfloor a \rfloor\)，\(z := 1 - 2(a - N)\)， \[ C_n(z) := \frac{1}{2^{2n}} \sum_{k=0}^{\lfloor 3n/4 \rfloor} \frac{d^{(n)}_k}{\pi^{2n - 2k} (3n - 4k)!} F^{(3n-4k)} (z), \] 再其中有 \[ F(z):=\frac{\cos \frac{\pi}{2}\left(z^{2}+\frac{3}{4}\right)}{\cos \pi z}, \] 以及 \[ \begin{aligned} d_{k}^{(n+1)}& =(3 n+1-4 k)(3 n+2-4 k) d_{k}^{(n)}+d_{k-1}^{(n)} & & n \geq 0 \text{ and }0 \leq k<\frac{3(n+1)}{4},\\ d_{k}^{(n)}& =0 && k<0 \text { or } k>\frac{3 n}{4},\\ d_{3 n}^{(4 n)}& =\lambda_{n} && n \geq 0, \end{aligned} \] 以及 \[ \begin{aligned} \lambda_{0} &=1, \\ \lambda_{n+1} &= \frac{1}{n+1}\sum_{k=0}^{n} 2^{4 k+1}\left|E_{2 k+2}\right| \lambda_{n-k} & (n \geq 0). \end{aligned} \] 以及 \(E_{k}\) 为 Euler 数。

（此处给 MathPix 打个广告，五毛不谢。）

当然这里的 \(n\) 取到无穷大上去了，所以我们需要截断。截断误差也由 [Gabcke] 给出：

[Theorem 3.2.2, Gabcke] 定义 \[ \begin{equation} R_K(t) := \frac{1}{\sqrt{a}} \sum_{K+1}^{\infty} \frac{C_n(z)}{a^n}, \label{eq:gabcke-error-term} \end{equation} \] 当 \(t \geq 200\) 时，有 \[ \begin{array}{ll} \left|R_{0}(t)\right|<0.127 t^{-3 / 4}, & \left|R_{5}(t)\right|<0.061 t^{-13 / 4}, \\ \left|R_{1}(t)\right|<0.053 t^{-5 / 4}, & \left|R_{6}(t)\right|<0.661 t^{-15 / 4}, \\ \left|R_{2}(t)\right|<0.011 t^{-7 / 4}, & \left|R_{7}(t)\right|<9.2 t^{-17 / 4}, \\ \left|R_{3}(t)\right|<0.031 t^{-9 / 4}, & \left|R_{8}(t)\right|<130 t^{-19 / 4}, \\ \left|R_{4}(t)\right|<0.017 t^{-11 / 4}, & \left|R_{9}(t)\right|<1837 t^{-21 / 4}. \\ \end{array} \]

这货真是太复杂了，看得我脑袋疼……原文还是德语，我就只能看懂公式……（尽管如此，还是比 David Platt 的 paper 容易理解。）

（劝退结束）

可以注意到截断误差 \(R_K(t)\) 对于小的 \(t\) 不太友好，所以对于小的 \(t\) 还得靠 Euler-Maclaurin……[LO] 里也注意到了这里，但是问题不大，因为要求 \(\pi(x)\) 的话，对于 \(t < x^{1/4}\) 的 \(t\) 我们只会 evaluate \(\tO(n^{1/4})\) 次，所以对于这些 \(t\) 用 Euler-Maclaurin ，复杂度也是 \(\tO(n^{1/2})\)，不影响整体复杂度。

[Gacbke] 的算法复杂度显然是 \(\tO(\sqrt{t})\) 的，被第一项给 dominate 了。之后我们会讨论如何优化掉这一项。这里我们只要注意到 \(\cos(\theta(t) - t \log n) n^{-1/2} = \Re(e^{i \theta(t)} n^{-1/2-it})\) ，所以这一项本质上还是在求 \(\sum\limits_{n=1}^N n^{-s}\) 。

说到具体实现的话。这个定理可以说非常简单明了清晰易懂了，可能最大的问题是要求一个函数的高阶导数。一个方法是手动展开，反正只需要常数项，另一个方法是手写符号计算，大炮打蚊子，还可以用维护幂级数的加减乘除以及 apply 上一个函数（反正只需要 apply 上一个 cos 就行了），应该是常数最小的方法之一了。

这里我们主要是求 \(Z(t)\) ，因为 critical line 上基本上只关心零点，所以 \(Z(t)\) 就足够了。当然你要是真的想求 \(\zeta(s)\) 也是可以的，[de Reyna] 做了和 [Gabcke] 类似的分析，连公式都是差不多的，我就懒得再抄一遍了，反正也没人看 →_→ 有人还实现了 de Reyna 的算法。另一种方法是由 \(\zeta(s) = e^{-i\theta(t)} Z(t)\)，所以我们还需要 \(\theta(t)\) 的逼近。有了 Stirling 公式，\(\theta(t)\) 的近似非常简单：

[Eq 2, Brent] \[ \begin{equation} \theta(t) = \frac{t}{2} \log \frac{t}{ 2 \pi e} - \frac{\pi}{8} + \sum_{j=1}^\infty \frac{(1 - 2^{1 - 2j})|B_{2j}|}{4j(2j-1)t^{2j-1}}. \label{eq:RS-zeta} \end{equation} \]

[Galway] 的 \(\zeta(s)\) 算法

TL;DR 那你就别读了 →_→

[Galway] 的 \(\zeta(s)\) 算法比较 general ，可以对不在 critical line 上的 \(s\) 也求值。

……我懒得写 argument 了，直接上结果：

[Theorem 7.3, Galway] 令 \(H(z) := (1 - e^{2 i \pi z})^{-1}\)，\(z_1 := N + \frac{1}{2} - \alpha e^{i\pi/4}\) 为 \(N \swarrow N+1\) 上任意一点，\(h\) 为任意复数满足 \(\frac{h}{|h|} = -e^{i \pi /4}\)，对于任意 \(0 < N_L \leq R < N_R\)，有 \[ \begin{equation} \mathcal{R}(s) = \sum_{n=1}^N n^{-s} + h \sum_{m=0}^M f(z_1 + mh; s) - \sum_{n=N_L}^N n^{-s} H\left( \frac{n - z_1}{h} \right) + \sum_{n = N+1}^{N_R} n^{-s} H\left( \frac{z_1 - n}{h} \right) + \mathcal{E}_\Sigma + \mathcal{E}_f, \label{eq:galway-R} \end{equation} \] 其中 \(\mathcal{E}_\Sigma\) 和 \(\mathcal{E}_f\) 为余项，依赖于 \(h, N, N_L, N_R, M\) 的取值。

后面的分析给了如下算法：

[Section 7.4, Galway]

  1. 计算 \(z_s := \sqrt{s / (2 \pi i)}\)，取 \(N = \lfloor \Re z_s - \Im z_s \rfloor\)；

  2. 解出最大的 \(\alpha_1\)，使得 \(|f(z_1; s)| \le \epsilon\) 且 \(\alpha_1 > 0\)，取 \(z_1 = N + \frac{1}{2} + \alpha_1 e^{i \pi /4}\)；

  3. 解出最大的 \(\alpha_2\)，使得 \(|f(z_2; s)| \le \epsilon\) 且 \(\alpha_2 > 0\)，取 \(z_2 = N + \frac{1}{2} - \alpha_2 e^{i \pi /4}\)；

  4. 枚举 \(M\)，使得：

    1. 取 \(h = (z_2 - z_1) / M\)；

    2. 计算 \(z_L := \sqrt{(2h)^{-2} + s/(2\pi i)} + (2h)^{-1}\)，取 \(N_L = \lceil \Re z_L - \Im z_L\rceil\)；

    3. 计算 \(z_R := \sqrt{(2h)^{-2} + s/(2\pi i)} - (2h)^{-1}\)，取 \(N_R = \lfloor \Re z_R - \Im z_R \rfloor\)；

    4. 计算 \(\mathcal{E}_f\) 的上界：令 \(g(z; s) := z^{-s} e^{i\pi z^2}\)，有 \[ \left|\mathcal{E}_{f}\right| \ll\left|g\left(z_{\mathcal{L}} ; s\right) e^{-2 i \pi\left(z_{\mathcal{L}}-z_{1}\right) / h}\right|+\left|g\left(z_{\mathcal{R}} ; s\right) e^{2 i \pi\left(z_{\mathcal{R}}-z_{1}\right) / h}\right| . \]

    5. 找最小的正整数 \(M\) 使得 \(|\mathcal{E}_f| \leq \epsilon\)。

  5. 用以上过程找到的 \(N, z_1, M, h, N_L, N_R\) 来计算 \(\eqref{eq:galway-R}\)。

[Galway] 没有给出一个 explicit error bound。他只在 [Section 7.7, Galway] 中，给了一个 asympototic error bound：

[Section 7.7, Galway] \[ \begin{aligned} |\mathcal{E}_f| & = O(|f(z_1; s)| + |f(z_2; s)|),\\ |\mathcal{E}_\Sigma| & = O(\left|g\left(z_{\mathcal{L}} ; s\right) e^{-2 i \pi\left(z_{\mathcal{L}}-z_{1}\right) / h}\right|+\left|g\left(z_{\mathcal{R}} ; s\right) e^{2 i \pi\left(z_{\mathcal{R}}-z_{1}\right) / h}\right| ). \end{aligned} \]

所以我也给不出分析 233333

在第 2 步和第三步涉及到最大化一个目标函数。在我实验中，这个函数几乎是单调的，所以二分一下就行了。

复杂度的话，明显可以看到是 \(\tO(N + M + N_R - N_L)\)。显然 \(N = \tO(t^{1/2})\)。[Section 7.6, Galway] 说要想做到 \(\epsilon\) 近似的话， \(M = \tO(\ln^{3/2} \epsilon)\)，而实验证明 \(N_R - N_L\) 被 \(M\) dominate 了，所以复杂度就是 \(\tO(t^{1/2})\)，而且瓶颈也在 \(\eqref{eq:galway-R}\) 第一项，和 Riemann-Siegel 公式一样。在下一大节里，我们讲讨论如何优化 \(\sum\limits_{n=1}^N n^{-s}\)。

[Zetafast]

这是我在网上搜到的一篇 paper，方法似乎和 Riemann-Siegel 公式不同，但是我没仔细看了。另外我始终觉得他的复杂度估的不对：\(D(N, s)\) 里面需要求 \(\lambda v N\) 次 \(Q(v, \cdot)\)，如果暴力求 \(Q(v, \cdot)\) 的话，总复杂度就变成 \(O(\lambda v^2 N) = O(v \sqrt{v \tau})\) 而不是 claim 的 \(O(\sqrt{vt})\)。作者 Kurt Fischer 我也搜不到人，但是 arxiv 提供的源代码里，这人自称是日本津山工业高等专门学校的。而这 paper 质量也还行，不像是出自民科之手。这里有 pari/gp 的实现。

\(\zeta(s)\) 预处理快速求值

到这里，用 Riemann-Siegel 公式，即使没有预处理，我们也可以得到一个 \(\tO(x^{3/5})\) 的算法来求 \(\pi(x)\)：只要好好 balance 一下暴力和积分的复杂度就行了。例如我们用 [Galway] 的思路，取 \(\lambda = \tO(x^{-2/5})\)，这样积分上限 \(T =\tO(x^{2/5})\)，每次用 \(\tO(x^{-1/5})\) 的时间求一次 \(\zeta(s)\) ，整体复杂度就是 \(\tO(^{3/5})\)。但是作为 bibi 选手，我们怎么能停在这里呢！反正都是理性愉悦，不如将嘴巴精神贯彻到底！

前文提到的 [Gabcke] 和 [Galway] 都是基于 Riemann-Siegel 公式，复杂度为 \(\tO(t^{1/2})\)，且瓶颈都是在求 \(F(t; N) :=\sum\limits_{n=1}^N n^{-(\sigma + it)}\) 这一项，这里的 \(\sigma\) 是固定的。瓶颈之外，两者算法都是 \(\text{poly}\log \epsilon^{-1}\) 级的。这里我们就来介绍一下如何通过预处理来快速求 \(F(t; N)\)。为了表达简单，接下来我就省略 \(F(t; N)\) 里的 \(N\) 了。

我 survey 到的算法均和 [OS] 类似，里面的 O 是之前提到的 A.M. Odlyzko，S 是 Schönhage，他提出的 Schönhage-Strassen 算法已经被 GMP 实现出来了。这类算法的思路为：我们希望在 \(\tO(N)\) 的时间内算出 \(F(t_0), F(t_0 + \delta), \dots, F(t_0 + M \delta)\)，其中 \(M, \delta\) 是两个我们可以控制的参数，然后再设计不同的算法用预处理的 \(F\) 快速计算满足 \(t_0 \leq t \leq t_0 + M \delta\) 的任何一个 \(F(t)\)。所以我们会有两部分：一部分如何预处理，一部分如何在线回答。

[FKBJ] 的预处理

TD; DR: 找到一个次数合适的多项式 \(P(t)\) 来逼近 \(e^{it}\) ，使得 \(\eqref{eq:poly-approx-err}\) 不超过 \(\epsilon\)。然后用 \(\eqref{eq:fft-f}\) 计算 \(f(x; k)\) ，再用快速 Fourier 变换计算 \(\hat f(x; k)\)，最后通过 \(\eqref{eq:final-G}\) 来计算所有的 \(G(t)\)。

他们号称是对 [Odlyzko] 的改进。我们考虑一个更 general 的问题：给定复数列 \(a_{1, \dots, N}\) 和实数列 \(\gamma_{1, \dots, N}\)，我们希望求对 \([-T, T]\) 里的所有整数 \(t\) 同时求 \[ G(t) := \sum_{n=1}^N a_n e^{i \gamma_n t}. \] 取 \(a_n = n^{-(\sigma + i t_0)}, \gamma_n = -\log n\)，这时的 \(G(t)\) 就是我们的目标 \(F(t)\) 了。

首先，我们令 \(R = 2^r > T\) 为最小的大于 \(T\) 的 2 的幂，于是对于任意 \(\gamma_n\)，都都可以找到一个 \(\delta_n = \gamma_n - \frac{2 \pi d_n}{R}\) 且 \(|\delta_n| \leq \frac{\pi}{R}\)，\(d_n\) 为整数。然后，我们可以找到一个 \(K\) 次多项式 \(P(t) = \sum\limits_{k=0}^K c_k t^k\) 来在 \(\left[\frac{-\pi T}{R}, \frac{\pi T}{R} \right]\) 上逼近 \(e^{it}\)。

于是我们有 \[ G(t) = \sum_{n=1}^N a_n e^{i \gamma_n t} = \sum_{n=1}^N a_n e^{2\pi i d_n t/R + i\delta_n t} = \sum_{n=1}^N a_n e^{2 \pi i d_n t / R} P(\delta_n t) + \mathcal{E}, \] 其中 \(\mathcal{E}\) 为用 \(P(t)\) 近似 \(e^{it}\) 带来的误差。

继续展开 \(P(\delta_n t)\) ： \[ G(t) = \sum_{k=0}^K c_k t^k \sum_{n=1}^N (a_n \delta_n^k) e^{2 \pi i d_n t/R} + \mathcal{E}. \] 我们惊奇的发现，里面这个 \(\sum\limits_{n=1}^N\) 长得和 Fourier 变换好像……不妨令 \[ \begin{equation} f(x; k) := \sum_{n \in [N]: d_n \equiv x \bmod{R}} a_n \delta^k_n, \qquad \hat f(x; k) := \sum_{x=0}^{R-1} f(x; k) e^{2 \pi i x/R}, \label{eq:fft-f} \end{equation} \] 那不就 \(\hat f\) 是 \(f\) 的离散 Fourier 变换嘛……

于是有 \[ \begin{equation} G(t) = \sum_{k=0}^K c_k t^k \hat f(x; k) + \mathcal{E}. \label{eq:final-G} \end{equation} \] 误差分析的话：这个比较简单，我上我也行， \[ \begin{equation} |\mathcal{E}| \leq \mathcal{E}_{e} = \sum_{i=1}^N |a_n| \max_{|t| \leq \pi T/R} |e^{it} - P(t)|, \label{eq:poly-approx-err} \end{equation} \] 所以只要 \(P(t)\) 和 \(e^{it}\) 足够靠近，问题就不大。用 \(e^{it}\) 的泰勒展开，只要 \(O(\log \epsilon^{-1})\) 项就能逼近到 \(|\mathcal{E}| \leq \epsilon\)。

不难看出，预处理复杂度是 \(\tO(N)\) 的。

[OS] 的在线计算

懒得看了……有兴趣的朋友也可以参考 [Section 2.3, Gourdon]。我在这里只是道听途说的参考 [Odlyzko]。

[OS] 的思路是：我仍然预处理出一个格点上的 \(F(t; N)\)。不仅如此，我还预处理出格点上的 \(F(t; N)\) 上的高阶导数： \[ F^{(m)}(t; N) = \sum_{n=1}^N n^{-\sigma-it} (-i \ln n)^m = \sum_{n=1}^N (n^{-\sigma} (\ln n)^m) e^{-i t \ln n}, \] 后者仍然是 \(G(t)\) 的形式，所以可以按照之前的方法一起处理。之后只要格点足够密，想求任意一个 \(F(t; N)\)，我们就可以在最近的一个格点处做泰勒展开。

[Gourdon] 和 [OS] 的在线计算

TL;DR: 确定 \(\beta, c, k_0, M\) 等参数，计算格点 \(\{t_0, t_0 + \pi / \beta, \dots, t_0 + M \pi / \beta\}\) 上的所有的 \(F(t; N)\)，然后用 \(\eqref{eq:gourdon}\) 计算 \(\eqref{eq:G-def}\) 定义的 \(G\)，最后用 \(\eqref{eq:decomp-F}\) 计算 \(F(t; N)\)。

现在我们假设我们能快速计算任何格点上的 \(F(t; N)\)，现在我们考虑如何用格点上的 \(F(t; N)\) 计算非格点上 \(F(t; N)\)。

[Odlyzko] 和 [Gourdon] 的方法差不多。相对而言，我更喜欢 [Gourdon] 的表述。由于 [Gourdon] 这块写的非常精妙，我这差不多就是他的翻译。（注：这里小节里请忘记上个小节的 \(G\) 函数，那是个临时变量。）

首先我们有如下定理：

[Whittaker–Shannon interpolation formula] 若一个函数 \(G(x)\) 可以表达为 \[ G(t) = \int_{-\tau}^\tau g(x) e^{i xt} \d x, \] ，定义 \(\sinc(x) = \frac{\sin(x)}{x}\)，则有 \[ G(t) = \sum_{n \in \mathbb{Z}} G\left( \frac{n\pi}{\tau} \right) \sinc(\tau t - n\pi). \]

可以看到，只要给定格点上的 \(G\) 值，Whittaker–Shannon 公式就能给出插值出整个函数。但是这样有一个问题：收敛太慢了。现在的差值公式只用 \(G\left(\frac{n\pi}{\tau} \right)\) 里的点，我们是否能通过增加格点的密度来加速收敛呢？

于是 [Gourdon] 给出了如下性质：

[Proposition 3, Gourdon] 给定两个复函数 \(G\) 和 \(h\) ，如果对于任意 \(z\) 均有 \(|G(z)| = O(e^{\tau |\Im z|})\)，\(|h(z)| = o(e^{(\beta - \tau)|\Im z|})\) 且 \(h(0) = 1\)，其中 \(\beta > \tau > 0\) 为两常数，则有 \[ \begin{equation} G(t) = \sum_{n \in \mathbb{Z}} G\left( \frac{n \pi}{\beta} \right) h(x - n \pi / \beta) \sinc(\beta x - n \pi). \label{eq:gourdon} \end{equation} \]

证明也非常简单。我都忍不住要抄一发……

[Proof of Proposition 3, Gourdon] 固定 \(x\)，考虑以下路径积分 \[ \int_C G(z) \frac{\beta}{\sin \beta z} \frac{h(x - z)}{x - z} \d z, \] 其中 \(C\) 是 0 为圆心，半径为 \((n+1/2) \pi / \beta\) 的圆。不难看出，随着 \(n\) 的增大，\(\frac{G(z) h(x - z)}{\sin \beta z} = \frac{o(e^{\beta |\Im z|})}{\Theta(e^{\beta |\Im z|})}\) ，故当 \(n \to \infty\) 时，该积分值为 0。而在我们取 \(n \to \infty\) 的极限的过程中，所有极点都被包含了进去。明显可以看出，这个积分的极点为 \(x\) 以及所有整数 \(n\) 对应的 \(n \pi / \beta\)。由留数定理得 \[ \text{Res}(x) + \sum_{n \in \mathbb{Z}} \text{Res}(n\pi/\beta) = 0, \] 计算可得 \[ \text{Res}(x) = -G(z) \frac{\beta}{\sin \beta z}, \quad \text{Res}(n \pi / \beta) = G\left( \frac{n \pi }{\beta} \right) \frac{1}{\cos (n\pi)} \frac{h(x - n \pi / \beta)}{x - n \pi / \beta}, \] 代入整理即可。

总之，有了 \(\eqref{eq:gourdon}\)，我们就可以设计一个好的 \(h\)，通过增加格点密度来加速收敛。关于 \(h\) 的设计，[Gourdon] 和 [Odlyzko] 不同。[Gourdon] 用的是 \(h(t; \beta, \tau, m) = \sinc^m\left( (\beta - \tau)t/m \right)\)，而 [Odlyzko] 则使用了 \[ h(t; \beta, \tau, c) = \frac{c}{\sinh c} \frac{\sinh \sqrt{c^2 - \eps^2 t^2}}{\sqrt{c^2 - \eps^2t^2}}, \quad \eps:= (\beta - \tau) / 2. \] 最后就是这个 \(G\) 到底是啥了。[Gourdon] 和 [Odlyzko] 都说：函数 \[ \begin{equation} G(t; k_0, k_1) = e^{-i \alpha t} \sum_{n=k_0}^{k_1} n^{-1/2 -it}, \quad \alpha := \ln \sqrt{k_0 k_1}, \label{eq:G-def} \end{equation} \] 对应的 \(\tau\) 为 \(\ln k_1 - \ln k_0\)，他们说是就是吧。这里我加了求和上下限，因为这样可以减小 \(\tau\)，减小格点密度，反正每次就暴力 \(k_0 - 1\) 次。 \[ \begin{equation} F(t; N) := \sum_{n=1}^N n^{-1/2-it} = \sum_{n=1}^{k_0 - 1} n^{-1/2 - it} + G(t; k_0, N). \label{eq:decomp-F} \end{equation} \] 当然对于 \(\eqref{eq:gourdon}\) 我们也不可能把所有的 \(n\) 都用上，必须截断。这样就会有截断误差，然而我懒得分析截断误差了……而且我们求格点上的 \(G(t; k_0, k_1)\) 也有误差需要考虑。[Gourdon] 上的分析看得我云里雾里，我就不总结了。如果用了 [Odlyzko] 的 \(h\)，那么我们只需要 \((\beta - \tau) t / 2 \leq c\) 的 \(t\) 就行了（因为我也不知道 \(\eps t > c\) 会发生什么……变成虚数？），而且 \(\sinh c \approx e^c / 2\)，\(h(t; \beta, \tau, c)\) 收敛的挺快的。但是这个算法里面有一堆参数，例如 \(\beta, c, k_0\)，我也不知道怎么选。[Equation 4.4.21, Odlyzko] 给出了他用的参数，可以当个参考。

我看到 [Odlyzko] 的这个 \(h\) 很眼熟呀，长得和 [FKBJ] 里的 kernel 好像。这个 kernel 起源是 [Logan]，我后来一时兴起就去看了一下 [Logan] 的 paper ，发现 Abstract 里面有这么一句话

...This result confirms a conjecture of J.C. Lagarias and A.M. Odlyzko ...

这不就是 [LO] 吗……难怪 [FKBJ] 会用这个 kernel 2333333

[Platt2] 的算法

没看懂……看起来他不是基于 Riemann-Siegel 公式的。就先放着吧……

下期预告

好了到这里，我们就可以均摊 \(\tO(1)\) 地对一个固定的 \(\sigma\)求 \(\zeta(\sigma + it)\) 了，用最开始的 [LO] 算法或者 [Galway] 里的算法，我们就可以做到 \(\tO(x^{1/2})\) 的时间内求 \(\pi(x)\) 了，算是说完了一个故事了，赶紧棺材板盖上完成从入门到入土。但是 [FKBJ] 和 [Platt] 的故事并没有完全结束，因为他们依赖于找到 \(\zeta(s)\) 的非平凡零点。虽然我们现在已经学会如何求 \(\zeta(s)\) 了，但是怎么找零点显然还需要一些额外的知识，包括：如何找到包含零点的区间？如何确定零点已经被找完了，没有遗漏？预知后事如何，请听下回分解……（我真写的要吐了……）

References

  • [LO] : Lagarias and Odlyzko, “Computing \(\pi(x)\): an analytic method.”
  • [OS] : Odlyzko and Schönhage, "Fast algorithms for multiple evaluations of the Riemann zeta function"
  • [Galway] : William Floyd Galway, "Analytic Computation of the Prime-Counting Function",
  • [FKBJ] : Franke et al., “A Practical Analytic Method for Calculating \(\pi (x)\).”
  • [Platt] : David Platt, “Computing \(\pi(x)\) Analytically.”
  • [Platt2] : David Platt, “Isolating some non-trivial zeros of zeta”
  • [de Reyna] : Juan Arias de Reyna, "High Precision Computation of Riemann's Zeta Function by the Riemann-Siegel Formula, I"
  • [Gabcke] : Wolfgang Gabcke, "Neue Herleitung und explizite Restabschätzung der Riemann-Siegel-Formel"
  • [Zetafast] : Kurt Fischer, "The Zetafast algorithm for computing zeta functions"
  • [Brent] : Richard P. Brent, "On asymptotic approximations to the log-Gamma and Riemann-Siegel theta functions"
  • [Gourdon] : Xavier Gourdon, "The \(10^{13}\) first zeros of the Riemann Zeta function, and zeros computation at very large height"
  • [Logan] : Benjamin F. Logan, "Bounds for the tails of sharp-cutoff filter kernels"