$\newcommand{\tO}{\widetilde O} \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”