\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.”