0%

解析方法数质数:从入门到入土(下)

\(\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"