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

Posted on In Disqus: Word count in article: 5.2k Reading time ≈ 19 mins.

\(\newcommand{\eps}{\epsilon} \newcommand{\tO}{\widetilde O} \newcommand{\d}{\mathrm d} \newcommand{\hphi}{\widehat \phi} \newcommand{\hPhi}{\widehat \Phi}\) 前一段时间某群里有人发了一个链接，说的是有人又整出一个新的算 \(M(x) = \sum\limits_{n=1}^x \mu(n)\) 的算法了，时间复杂度 \(\tO(x^{3/5})\)，空间复杂度 \(\tO(x^{3/10})\)。然后有人提到目前最快的算 \(M(x)\) 的算法是用解析数论的方法，和 \(\pi(x)\) 类似，时间复杂度 \(\tO(x^{1/2})\)，空间复杂度 \(\tO(x^{1/4})\) 。那几篇论文都是针对 \(\pi(x)\)，然后说可以拓展到 \(M(x)\)，我一时兴起就去学习了一下，整一点理性愉悦的东西。

在这个过程中，我读了好几篇 paper，有几篇没读懂，这里差不多就是对那几篇 paper 的 survey 吧……由于这个领域属于 analytic number theory，而我不懂复分析，所以……几乎所有核心定理的证明我都看不懂……所以由于个人能力原因，文章中可能存在不严谨的地方或者错误，大家意思意思一下吧 →_→

由于内容比较多，所以我可能会拆成好几篇文章来写（又可以水了呢 →_→）。

历史

这里我首先介绍一下解析方法数质数这个小领域的历史。首先是 1987 年 [LO] 整了一个算法，说是可以在 \(\tO(n^{1/2})\) 的时间内算 \(\pi(x)\)，但是由于常数过大，他们并没有实现出来……后来 William Galway 在 UIUC 数学系读博的时候一直研究 [LO] 的算法，做了点小改进，补充了里面的技术细节。然而，Galway 自己也没实现出来……再后来 [FKBJ] 用另外一种思路做了点小改进，然后真的实现出来了！[FKBJ] 在假定 Riemann hypothesis 的前提下，花了几千个 CPU hours 算出来 \(\pi(10^{24})\)。后来 [Platt] 在 [Galway] 的基础上也做了点小改进，然后实现了代码！他也算了 \(\pi(10^{24})\)，在没有 assume RH 的前提下确认了 [FKBJ] 的结果。再后来 [Büthe] 在 [FKBJ] 的基础上又做了点小改进，算出了 \(\pi(10^{25})\)。这个 Büthe 也是 FKBJ 里的 B。再后来的事我就没注意了……

主要思路

这里我们先约定：所有字母 \(p\) 均表示素数。于是 \(\pi(x)\) 可以表示为 \(\pi(x) = \sum\limits_{p \leq x} 1\)，尽管这样做什么用也没有……文章中所有的 \(i\) 均为 \(\sqrt{-1}\)。复杂度里的 \(\widetilde O\) 和 \(\widetilde \Omega\) 里均省略一个 \(\text{poly}\log n\)，因为我实在懒得算有多少个 \(\log\) 了……文中的 \(\log\) 和 \(\ln\) 混用了。

这一切要从 [LO] 说起……作者之一的 Andrew M. Odlyzko 接下来会出现很多次。但是由于 [LO] 写的比较意识流，我沿用 [Galway] 的思路。

万恶之源是 Mellin 变换：给定一个函数 \(\phi(s)\)，它的 Mellin 变换为 \[ \hphi(s):= \int_0^\infty \phi(u) u^{s - 1} \d u. \] 然后就有定理 [Theorem 1.4, Galway] 告诉我们：对于任意一个定义域内的 \(\sigma\)，我们有 \[ \begin{equation} \frac{\phi(u+) + \phi(u-)}{2} = \frac{1}{2\pi i} \int_{\sigma - i \infty}^{\sigma + i\infty} \hphi(s) u^{-s} \d s. \label{eq:mellin-inv} \end{equation} \] 根据我的猜测，这里的 \(\phi(u+)\) 和 \(\phi(u-)\) 应该是指 \(\lim\limits_{\delta \to 0^+}\phi(u + \delta)\) 和 \(\lim\limits_{\delta \to 0^-} \phi(u + \delta)\)。如果 \(\phi(u)\) 性质没那么差，我们就有 \(\phi(u) = \frac{\phi(u+) + \phi(u-)}{2}\)。注意这个 \(\phi(u)\) 可以是不连续的。这样，\(\eqref{eq:mellin-inv}\) 可以看成是 Mellin 变换的逆变换。给定一个序列 \(\{a_n\}_{n \in \mathbb{N}^+}\)，我们对逆变换稍作处理： \[ \begin{equation} \sum_{n \geq 1} a_n \phi(n) = \frac{1}{2\pi i} \int_{\sigma - i \infty}^{\sigma + i\infty} \hphi(s) \sum_{n \geq 1} a_n n^{-s} \d s. \label{eq:weighted-mellin} \end{equation} \] 左边看起来平淡无奇，但是右边的 \(\sum\limits_{n \geq 1} a_n n^{-s}\) 不免让人遐想联翩：这不就是 Dirichlet 级数嘛！例如我们令 \(a_n = \mu(n)\)，我们有 \(\sum\limits_{n \geq 1} a_n n^{-s} = \frac{1}{\zeta(s)}\) ，那么两边就可以化简成 \[ \sum_{n \geq 1} \mu(n) \phi(n) = \frac{1}{2\pi i} \int_{\sigma - i \infty}^{\sigma + i \infty} \hphi(s) \zeta(s)^{-1} \d s. \] 如果我们令 \(a_n\) 表示 \(n\) 是否为素数，那么这个 Dirichlet 级数还是不好求。没关系，我们可以令 \(a_{p^m} = \frac{1}{m}\)，其它的 \(a_n = 0\) ，这样就有： \[ \sum_{n \geq 1} a_n n^{-s} = \sum_{p^m} \frac{1}{m} p^{-ms} = \sum_{p}\sum_{m \geq 1} \frac{1}{m} p^{-ms} = \sum_p \ln (1 - p^{-s})^{-1} = \ln \prod_{p} \frac{1}{1 - p^{-s}} = \ln \zeta(s). \] 还是可以继续的嘛……于是我们有： \[ \begin{equation} \sum_{p^m} \frac{1}{m} \phi(n) = \frac{1}{2\pi i} \int_{\sigma - i \infty}^{\sigma + i \infty} \hphi(s) \ln \zeta(s) \d s. \label{eq:pi-star-mellin} \end{equation} \] 但是这样我们就不再求 \(\pi(x)\) ，转而求 \[ \begin{equation} \pi^*(x) = \sum_{p^m \leq x} \frac{1}{m} = \sum_{m \geq 1} \frac{1}{m} \pi(\sqrt[m]{x}). \label{eq:pi-pi*} \end{equation} \] 注意到 \(m > 1\) 时，\(\sqrt[m]{x}\) 衰减的很快，可以直接暴力了，所以由 \(\pi^*(x)\) 推 \(\pi(x)\) 问题不大，时间复杂度 \(\tO(x^{1/2})\)。比较 \(\eqref{eq:pi-pi*}\) 和 \(\eqref{eq:pi-star-mellin}\)，自然而然，我们就想让 \(\phi(n) = [n \leq x]\) ，然而这让 \(\hphi(s)\) 难算，右边的积分积不出来……所以我们需要一个平滑一点的函数（\([n \leq x]\) 这个函数太阶梯了），满足 \(\phi(n)\) 在 \((-\infty, x_l]\) 非常接近于 1，\([x_r, \infty)\) 上非常接近于 0 ， \(x_l \leq x \leq x_r\) 且 \(x_r - x_l\) 不太大，这样我们还是可以强行算积分，然后暴力把 \([x_l, x_r]\) 里的不同减掉即可。详细说来，我们想找一个 \([n \leq x]\) 的近似 \(\phi(u)\)，然后计算 \[ \begin{equation} \pi^*(x) = \sum_{p^m} \frac{1}{m} [p^m \leq x] = \frac{1}{2 \pi i} \int_{\sigma - i \infty}^{\sigma + i \infty} \hphi(s) \ln \zeta(s) \d s + \sum_{p^m} \frac{1}{m} \left( [p^m \leq x] - \phi(p^m) \right). \label{eq:key} \end{equation} \] 右式分为两部分，第一部分为积分，第二部分为用 \(\phi(n)\) 来近似 \([n \leq x]\) 带来的误差。我们希望 \(\phi(n)\) 只在 \([x_l, x_r]\) 上有可能和 \([n \leq x]\) 有较大差别，在 \((-\infty, x_l]\) 和 \([x_r, \infty)\) 上和 \([n \leq x]\) 的差距小到可以忽略。请记住 \(\eqref{eq:key}\)，这个等式是这一系列算法的核心，所有的算法都围绕等式右边来设计，特别是里面的积分这一项。

到这里，我们总结一下这个算法的思路：

  1. 设计一个好的 \(\phi(u)\) 使得 \(\eqref{eq:key}\) 右边两项均可以快速计算；
  2. 用 \(\eqref{eq:key}\) 计算 \(\pi^*(x)\)，时间复杂度依赖于 \(\phi(u)\)；
  3. 然后用 \(\eqref{eq:pi-pi*}\) 求 \(\pi(x)\)，这一步时间复杂度 \(\tO(x^{1/2})\)。

接下来，[LO] 和 [Galway] 分别设计了自己的 \(\phi(u)\) 以及计算 \(\eqref{eq:key}\) 的算法，复杂度均为 \(\tO(x^{1/2})\)。[Platt] 用的是 [Galway] 的 \(\phi(u)\)，但是设计了新的算法来算 \(\eqref{eq:key}\) 的积分。

[LO] 的 \(\phi(u)\)

首先我们给出定义：固定 \(0 < w < x, k\) ，令 \(f(u; k) = c^{-1} u^k (1 - u)^k\) 其中 \(c = \int_0^1 u^k (1-u)^k \d u\) 为归一化常数，则 \(\phi(u)\) 定义为 \[ \phi(u; x, w, k) = \begin{cases} 1 & u \leq x - w, \\ 0 & u \geq x, \\ \int_0^{(x - u)/w} f(v; k) \d v & x - w < u < x. \end{cases} \] 其 Mellin 变换为 \[ \hphi(s; x, w, k) = s^{-1} \int_0^1 (x - vw) f(v; x, w, k) \d v. \] 怎么求 \(\phi(u; x, w, k)\) 我没看懂，paper 里面这一页的排版直接劝退……我猜是 Euler-Maclaurin 算个积分就行了，或者可以整出一个 closed-form 出来……所以只要 \(x - w \leq x \leq x\) ，算 \(\eqref{eq:key}\) 第二项就可以在 \(\tO(w)\) 内完成。现在问题的重点就是如何求 \(\eqref{eq:key}\) 里的积分。方便起见，我们定义 \(\Psi(s; x, w, k) = \hphi(s; x, w, k) \log \zeta(s)\)。再注意到 \(\hphi(\bar s; x, w, k) = \overline{\hphi(s; x, w, k)}\)，\(\zeta(\bar s) = \overline{\zeta(s)}\)，故有 \(\Psi(\bar s; x, w, k) = \overline{\Psi(s; x, w, k)}\)。由于积分的对称性，我们只需要考虑 \(\Re\Psi(s; x, w, k)\) 的积分即可。

我们把这个积分的上下限从 \((-\infty, \infty)\) 截断为 \([-T, T]\)： \[ \frac{1}{2 \pi i} \int_{\sigma - i \infty} ^{\sigma + i \infty} \Re \Psi(s; x, w, k) \d s \approx \frac{1}{2 \pi i} \int_{\sigma - i T} ^{\sigma + i T} \Re \Psi(s; x, w, k) \d s + \mathcal{E}_T, \] 其中 \(\mathcal{E}_T\) 为误差项。[Eq 3.12, LO] 给出了 \(\mathcal{E}_T\) 的 bound： \[ |\mathcal{E}_T| = \tO(x^{3/2+k} (Tw)^{-k}), \] 所以 \(T = \tO(x^{1+3/(2k)}/w)\) 就可以把截断误差控制住。但是即使截断了，这个积分仍然不好求。[LO] 给出的方法是 Euler-Maclaurin 求和公式：

[Lemma 3, LO] 如果 \(f(t)\) 在 \([0, T]\) 上有 \(2m\) 次导数，那么对于任意 \(n \in \mathbb{N}\) 有 \[ \begin{equation} \begin{aligned} \int_0^T f(t) \d t & = \frac{T}{n} \left[\sum_{j=0}^n f\left(\frac{Tj}{n}\right) - \frac{g(0) + g(T)}{2} \right] \\ & \hspace{4em} - \sum_{j=1}^{m-1} \frac{B_{2j}}{(2j)!} \left( \frac{T}{n} \right)^{2j-1} \left[g^{(2j-1)}(T) - g^{(2j-1)}(0) \right] + R_{2m}, \end{aligned} \label{eq:euler-maclaurin} \end{equation} \] 其中 \(B_{2i}\) 为 Bernoulli 数，\(R_{2m}\) 为误差项 \[ R_{2m} \leq n \frac{|B_{2m}|}{(2m)!} \left( \frac{T}{n} \right)^{2m+1} \max_{t \in [0, T]} \left| g^{(2m)}(t) \right|. \]

意会一下，我们想通过在 \([0, T]\) 内均匀撒 \(n\) 个点来求 \(\int_0^T f(t) \d t\)，[Lemma 3, LO] 给出了这样的做的误差。可以看到，\(n\) 越大，误差约小，误差项 \(R_{2m}\) 则和 \(f^{(2m)}(t)\) 的最大值相关。[Lemma 4, LO] 给出了 \(\Re\Psi(\sigma + it; x, w, k)\) 的 \(m\) 次导数的 bound：

[Lemma 4, LO] \[ \left| \frac{\d^m}{\d t^m} \Re \Psi(\sigma + it; x, w, k) \right| = O\left( m! ((\sigma-1)/2)^{-m} x^{(\sigma+2)/2+k} \right). \]

故我们取 \(\sigma = 3/2\)，\(m = k^2\)， \(n = \tO(x^{1+3/k}/w)\) 就能把误差控制在合理范围内。（这句话我瞎 bibi 的，没算。）

很明显， \(\eqref{eq:key}\) 第二项的复杂度只取决于 \(w\)，而假定 \(\Psi(s; x, w, k)\) 及其 \(m\) 阶导数能在 \(\tO(1)\) 的时间内（\(\tO\) 内省略了 \(m\)）求出来的话，第一项积分的复杂度为 \(\tO(n) = \tO(x^{1+3/k} / w)\)，取一个足够大的 \(k\) 就能使得这一项复杂度变为 \(\tO(n/w)\)。故为了平衡两边复杂度，我们取 \(w = \tO(n^{1/2})\) ，两边时间复杂度就都是 \(\tO(x^{1/2})\)。

最后，我们来总结一下 [LO] 算法是如何计算 \(\eqref{eq:key}\) 的：

  1. 第二项直接暴力，复杂度 \(\tO(w) = \tO(x^{1/2})\)；
  2. 第一项我们先把积分截断到 \(T = \tO(x/w) = \tO(x^{1/2})\) ，再用 Euler-Maclaurin 公式 \(\eqref{eq:euler-maclaurin}\) 采 \(n = \tO(x^{1+3/k}/w) = \tO(x^{1/2})\) 个点来计算。
  3. 整体时间复杂度 \(\tO(x^{1/2})\)。

下面是 paper 某一页的截图……诡异的排版真劝退。

[Galway] 的 \(\phi(u)\)

[Galway] 沿用了 [LO] 的思路，但是 [Galway] 换了一个 \(\phi(u)\)，他的 \(\phi(u)\) 定义如下： \[ \phi(u; x, \lambda) = \frac{1}{2} \text{erfc}\left(\frac{\ln (u / x)}{\sqrt{2} \lambda} \right). \] 对应的 Mellin 变换 \(\hphi(s)\) 为 \[ \hphi(s; x, \lambda) = e^{\lambda^2 s^2 / 2} \frac{x^s}{s}. \] 我们再来分析 \(\eqref{eq:key}\) 的第二项。注意到 \(\text{erfc}(x) \leq e^{-x^2}\)，所以 \(\phi(u; x, \lambda)\) 在 \(u\) 出了 \([x / (\lambda \sqrt{\log x}), x \lambda \sqrt{\log x}]\) 后衰减的非常快，所以这个函数几乎等价于阶梯函数了，所以我们就暴力把这个区间内的差求个和就行了，复杂度是 \(\tO(x\lambda)\)。

现在问题就是如何求 \(\eqref{eq:key}\) 里的积分了。和上文的分析一样，我们继续定义 \(\Psi(s;x, \lambda) = \hphi(s; x, \lambda) \ln \zeta(s)\)，于是 \(\eqref{eq:key}\) 中的积分可以表示为 \(\int_0^{\infty} \Phi(\sigma + it) \d t\)。[Galway] 也是用两步来近似这个积分，不过把截断和离散化换了个顺序： \[ \int_0^\infty \Re \Psi(\sigma + it; x, \lambda) \d t \approx h \sum_{k=0}^{\infty} \Re \Psi(\sigma + iht; x, \lambda) \approx h \sum_{k=0}^{T/h} \Re \Psi(\sigma + iht; x, \lambda). \] 第一步是把积分换成了求和，[Section 3.3, Galway] 证明了只要 \(h = O\left(\frac{1}{\ln (x/\varepsilon)}\right)\)，那么积分转求和这一步的误差就可以被控制到 \(\varepsilon\)。第二步是把无限求和转化成有限求和，我们有如下定理：

[Theorem 3.25, Galway] 若 \(x > 0, \lambda > 0, \sigma > 1, h > 0, T > 0\)，则有 \[ h \sum_{k > T/h} |\Re \Psi (\sigma + ikh; x, \lambda)| \leq \frac{1}{2} e^{\lambda^2 \sigma^2 / 2} \ln \zeta(\sigma) E_1(\lambda^2 T^2 / 2), \] 其中 \(E_1(z) = \int_z^\infty e^{-r} / r \d r\) 为指数积分。

所以只要 \(T = \widetilde\Omega\left(\lambda^{-1} \sqrt{\ln \epsilon^{-1}}\right)\) ，那么这一步的误差也能被控制到 \(\varepsilon\)。为了使得最后答案准确，我们只要使得这个积分误差不超过一个小常数即可，所以 \(\varepsilon\) 可以认为是常数（例如 0.1）。假定我们能 \(\tO(1)\) 时间内对 \(\Psi(s; x, \lambda)\) 进行一次求值，那么这个积分可以在 \(\tO(\lambda^{-1})\) 的时间内做到常数逼近。

综上，这个算法的总体时间为 \(\tO(x\lambda + \lambda^{-1})\)，所以令 \(\lambda = \widetilde \Theta(x^{1/2})\) 可以使得整体复杂度为 \(\tO(x^{1/2})\)。不过这里还埋了一个坑：如何在 \(\tO(1)\) 的时间内求一次 \(\Psi(s; x, \lambda)\)，特别是 \(\Psi\) 里面还涉及到 \(\zeta(s)\)。这事可不简单，可以让我再水一篇文章了 23333

[Galway] 还 argue 了一下某种意义上，他 propose 的这个 \(\phi(u)\) 是最优的。我没仔细看，以下是我的猜测：我们想要 \(\phi(u)\) 和阶梯函数比较像，才能快速算 \(\eqref{eq:key}\) 里的第二项；我们还想要 \(\hphi(s)\) 收敛的快，这样方便算第一项。如果我们把 \(\phi(u)\) 表达成 \(\phi(u) = \Phi(\ln(u/x) / \lambda)\) 这样的形式，其中 \(\Phi(\cdot)\) 可以是任意函数，那么 \(\hphi(s)\) 可以和 \(\Phi(\cdot)\) 的 Fourier 级数扯上关系（[Eq 2.37, Galway]），所以我们希望 \(\Phi\) 和 \(\Phi\) 的 Fourier 级数同时快速收敛。在这个意义下，他 propose 的这个 \(\phi(u)\) 是使得两者同时收敛最快的函数了。

[Platt] 的积分

首先我要吐槽一下 David Platt 的 writing，看得我实在太难受了，真的是 painful。他就给你列一堆 Lemma、Theorem，中间的联系都不写，每次我都得猜他这个 Lemma 想表达什么意思，是用来解决什么问题的。他写算法就给我罗列了一堆步骤，依次算 abcde 函数，看完了整个算法，我都不知道这个算法输出的是啥，和我们关心的东西有什么联系，abcde 各自有什么意义，为什么要这么做。[LO] 我能轻松跟上 high level idea，但是 David Platt 的 paper 我连 high level idea 都跟不上。要不是他 claim 自己真的码了代码实现了算法，我早就不看了……但也有可能是他文章就不是写给我这种不懂复分析的人看的……

废话不多说，鉴于我也看不懂 [Platt] 的思路，我就直接上公式了。

我直接引用里面最核心的一个定理 [Theorem 4.7, Platt]。

[Theorem 4.7, Platt] \[ \begin{equation} \begin{aligned} \frac{1}{2 \pi i} \int_{2 - i \infty}^{2 + i \infty} \hphi(s; x, \lambda) \ln \zeta(s) \d s & = \hPhi(1; x, \lambda) - \sum_{\rho} \Re \hPhi(\rho; x, \lambda) - \log 2\\ & \hspace{4em}+ \frac{1}{2\pi i} \int_{-1 - i \infty}^{-1 + i \infty} \hphi(s; x, \lambda) \ln (-\zeta(s)) \d s. \end{aligned} \label{eq:platt} \end{equation} \] 其中 \(\rho\) 为 \(\zeta(s)\) 的所有非平凡零点，且 \(\hPhi\) 定义为：对于 \(t > 0\) 有 \[ \hPhi(\sigma + i t; x, \lambda) = -\int_{\sigma + it}^{\sigma + i \infty} \hphi(s; x, \lambda) \d s, \] 对于 \(t < 0\) 有 \(\hPhi(\bar s; x, \lambda) = \overline{\hPhi(s; x, \lambda)}\)。

\(\hPhi\) 的准确定义你们可以自己去看 [Lemma 4.6, Platt]，我没看懂他的定义（再次吐槽 writing），瞎猜了一个，跑了跑代码似乎问题不大……

Anyway 假设我们的 \(\hPhi\) 定义没问题，我们就可以继续了。这个定理左边就是 \(\eqref{eq:key}\) 中的积分，只不过直接把 \(\sigma\) 变成了 2 。右边最后一个积分类似于 error term，[Lemma 4.5, Platt] 给出了 error term 的 bound：

[Lemma 4.5, Platt] \[ \left| \frac{1}{2 \pi i} \int_{-1 - i \infty}^{-1 + i \infty} \hphi(s; x, \lambda) \ln (-\zeta(s)) \right| \leq \frac{e^{\lambda^2 / 2} }{2 \pi x \lambda} \left( 5 \sqrt{2} \pi + 2 \lambda^{-1}\right). \]

所以只要 \(\lambda = \Theta((x \varepsilon)^{-1/2})\) ，这一项就能控制在 \(\varepsilon\)。

继续看 \(\eqref{eq:platt}\)，还有一件事需要处理：\(\zeta(s)\) 的非平凡零点有无数个，我们只能求一部分，怎么办？[Appendix A, Platt] 证明了，如果我们只考虑 \(|\Im \rho| \leq \tO(\lambda^{-1})\) 的 \(\rho\)，问题也不大。我实在懒得分析 [Lemma A.4, Platt] 里的 \(\alpha_{T_1}\) 和 \(N(T_1)\) 的关系了，这里就当这个差是 \(\tO(1)\) 了。这样截断误差可以被 bound 住。

另外还有一个问题：\(\hPhi(s; x, \lambda)\) 怎么求？[Section 6, Platt] 给出的方法很简单：泰勒展开。注意到，对于任意一个 \(\Re s_0 \neq 0\) 的点 \(s_0\)，先在 \(s_0\) 处化简（[Lemma 6.1, Platt]）： \[ \hphi(s_0 + ih; x, \lambda) = \hphi(s_0; x, \lambda) \exp\left( ih (s_0 \lambda^2 + \ln x) \right) \frac{\exp(-\lambda^2 h^2/2)}{1 + \frac{ih}{s_0}}. \] 注意到第一项和 \(h\) 无关，第二项可以认为是 \(\exp(c_1h)\) 其中 \(c_1 = i(s_0 \lambda^2 + \ln x)\)。不妨设 \(f(h; s_0, \lambda) = \exp(-\lambda^2 h^2 / 2) \left( 1 + \frac{ih}{s_0} \right)^{-1}\)，通过简单的分步积分，我们有 \[ \int \exp(c_1h) f(h; s_0, \lambda) \d h = \frac{\exp(c_1 h)}{c_1} \sum_{k=0}^{K-1} (-c_1)^{-k} f^{(k)}(h; s_0, \lambda) + (-c_1)^{-K} \int \exp(c_1 h) f^{(K)}(h; s_0, \lambda) \d h. \] 这样我们只要 bound 住 \(\max |f^{(K)}(h; s_0, \lambda)|\) 即可。 \(h\) 太大的话不太好 bound，所以我们希望 \(h\) 比较小。[Platt] 证明了，如果 \(\lambda h < 1\) 且 \(|h| < |s_0|\)，那么 \(\max |f^{(K)} (h; s_0, \lambda)|\) 是可以被 bound 的。所以我们做法是：把整个积分区间 \([t, \infty)\) 变成若干个不相交的小区间 \([t_0,t _1) \cup [t_1, t_2) \cup \cdots\)，每个区间内 \(h\) 最大是 \(t_j - t_{j-1}\)，每个区间内在 \(\sigma + i t_{j-1}\) 处展开，这样每个区间内都能保证 \(|\lambda h| < 1\) 和 \(|h| \leq |s_0|\) ，而且只需要 \(O(\log \lambda^{-1})\) 个区间就能求出 \([t, \lambda^{-1}]\) 的积分，再往上我们就直接舍去了……[Platt] 也没说这样做会引入多少误差……所以算一次 \(\hPhi(s; x, \lambda)\) 的时间可以控制在 \(\tO(1)\)。

（聪明的同学可以去看看 [Section 6, Platt]，看是我理解不行还是他写的不行 →_→ 不长，只有一页，notation 都是一样的。）

总结一下 [Platt] 的方法：

  1. [Platt] 提出了一个新的求 \(\eqref{eq:key}\) 中积分的方法，第二项暂不考虑，可以使用其余的算法；
  2. [Platt] 用 \(\eqref{eq:platt}\) 来求积分，最后一项误差项舍去；
  3. 在 \(\eqref{eq:platt}\) 中，我们需要计算 \(\hPhi\)，算法上文给出了；
  4. 未解决的问题：如何枚举出所有虚部不超过 \(T\) 的 \(\zeta(s)\) 的非平凡零点？而且为了控制误差，每个零点还得非常精确……

[FKBJ] 和 [Büthe]

看起来思路和 [LO] 不太一样……虽然还是看不懂，但是 [FKBJ] 上是数学上的看不懂，reference 了太多东西一下消化不了；David Platt 的是逻辑上的看不懂，上一段还是活蹦乱跳的一头猪，下一段怎么就变成了美味的小炒肉了呢？

Anyway 这两篇 paper 我就先鸽了，有时间再看 →_→

小结

这篇文章主要介绍了 [LO] 的思路及其设计的 \(\phi(u)\)，以及 [Galway] 和 [Platt] 对其思路的改进。中间还有很多坑没填，例如参数具体该怎么选，\(\tO\) 里到底省略了多少个 \(\log\) （手动 doge），以及最重要的：\(\zeta(s)\) 怎么求？无论是 [LO], [Galway] 还是 [Platt]，里面都需要疯狂 evaluate \(\zeta(s)\)，我都是假定 evaluation 可以在 \(\tO(1)\) 内搞定。但是具体怎么搞，这是一件非常不显然的事。而且，这几个算法都需要对 \(\zeta(s)\) 的 evaluation 有精度上的要求，下一篇文章，我就准备来说说，我 survey 到的几种求 \(\zeta(s)\) 的方法。

另外，虽然理论上时间复杂度 \(\tO(x^{1/2})\) 比组合算法（我知道最好的是这篇，复杂度 \(\tO(x^{2/3})\) ）优秀，但是实际跑起来嘛，怕不是被吊打（doge）……在 [Platt] 中有这么一段话：

[Platt] ...Assuming run time asymptotic to \(O(x^{2/3})\) and \(O(x^{1/2})\) for the combinatorial and analytic algorithms respectively, the crossover at which this implementation of the analytic algorithm would start to beat the combinatorial method would be in the region of \(x = 4 \cdot 10^{31}\).

但是现在我们也只算到了 \(\pi(10^{26})\)……

有兴趣的朋友可以看一下 StackExchange 上这个问题，DanaJ 的回答是比较全的了。

References

  • [LO] : Lagarias and Odlyzko, “Computing \(\pi(x)\): an analytic method.”
  • [Galway]: William Floyd Galway, "Analytic Computation of the Prime-Counting Function",
  • [FKBJ] : Franke et al., “A Practical Analytic Method for Calculating \(\pi (x)\).”
  • [Büthe] : Büthe, “An Improved Analytic Method for Calculating \(\pi(x)\).”
  • [Platt] : Platt, “Computing \(\pi(x)\) Analytically.”