超快筛法

· · 个人记录

大家都听过很多素数筛法，如时间复杂度为O(n\sqrt n)的朴素筛法，时间复杂度为O(n \log n \log n)的 Eratosthenes 筛法（埃拉托斯特尼筛法，又称埃氏筛），时间复杂度为O(n)的 Euler 筛法（欧拉筛，又称线性筛）[^\color{blue}1][sieves]。

但是，还有没有更快的筛法呢？

要解决传统的问题——求1\sim n的素数，筛法想低于线性非常困难。因为这就意味着这个筛法需要以一种类似函数的表达方法来存储1\sim n。而这种表达方法现在数学家们还没有找到。（通用的）

那如果我们求\pi(n)呢？

那么，神仙算法就有了。 为了各种方便，请先阅读 约定 。

因为the Meissel, Lehmer, Lagarias, Miller, Odlyzko method太难打了，在本文就记为A算法。（注意，是1985年那个版本） 在1988年M. DELEGLISE 和 J. RIVAT改进的算法记为B算法。

在本文里\mathbb Z表示正素数集（好像也没有人用负的），即\{x:x\in \mathbb N\wedge \forall\; 1 < y < x,y\nmid x\} 本文用\pi(x)表示不大于x的素数，即\left|[1,n]\cap\mathbb Z\right|。

如果未标明数的范围，默认是在正整数范围内的。

斜体表示注释。

  1. 什么是B算法

B算法，是一种快速求\pi (n)的算法。限于篇幅，本文将不会讲解此算法的原型——[A算法][final]。A算法在1985年被J. C. Lagarias， V. S. Miller 和 A. M. Odlyzko改进后的Meissel-Lehmer method，用于计算\pi(4\times 10^{16})。而在1996年M. Deleglise 和 J. Rivat将其优化，用于计算\pi (10^{18})。

B算法的时间复杂度是O\left(\dfrac{x^{2/3}}{\log^2x}\right)而空间复杂度为O\left(x^{1/3}\log^3x\log \log x\right)。[^\color{blue}2][paper]

  1. B算法是怎么运作的

本节仅为方便大家阅读，当然如果大家英语不错，看[原论文][paper]更好。

1.0 定义

$\phi(x,a)$表示$\leqslant x$的正整数里所有质因子都$>p_a$的数的个数。即 $\phi(x,a) = \left|\left\{n\leqslant x : n \in\mathbb N_+\land \forall p_j \mid n, p_j > p_a\right\}\right| 1.1 开始简单推导 容易得到一个公式（证明相信大家都会）： 可是$P_j(x,a)$的非零项并不是无穷多个的，因为 _（原文说的是$P_j(x,a)=0$，但是个人认为如果说是$a-1$会更好一点，因为后文要用）_ 然后我们定义一个正整数$y$，使得$x^{1/3}\leqslant y\leqslant x^{1/2}$。然后令$a=\pi(y)$。 然后， 由$(1.1.1)$得： 证明： 首先根据$\pi(x)$的性质，$p_{\pi(x)}\leqslant xx \therefore P_j(x,a)=0

然后我们得到

\begin{aligned} \phi(x, a) & = \sum^\infty_{j=0} P_j(x,a) \ & = \sum^2_{j=0} P_j(x,a) + \sum^\infty_{j=3} P_j(x,a) \ & = \sum^2_{j=0} P_j(x,a) \ & = P_0(x,a)+P_1(x,a)+P_2(x,a) \ & = 1+\pi(x)-a+P_2(x,a) \end{aligned} \therefore \phi(x, a)+a-1-P_2(x,a) = \pi(x) \therefore \pi(x) = \phi(x, a)+a-1-P_2(x,a)

所以，我们就把求\pi(x)的问题转化为求\phi(x,a)和P_2(x,a)的问题了。

1.2 求P_2(x,a)

那么，$p_a=p_{\pi(y)} \leqslant y <p\leqslant q$。 $\because p\in[y+1,\sqrt x],\quad q\in[p,x/p] 当$p\in[y+1,\sqrt x]$时，$\dfrac x p\in\left[1,\dfrac x y\right]$。 $\therefore$我们可以用埃氏筛筛了$\left[1,\dfrac x y\right]$，然后对于每个$\left[y+1,\sqrt x\right]$的素数$p$，求$\pi\left(\dfrac x p\right)-\pi(p)+1$的值并求和就能算出$P_2(x,a)$。时间复杂度为$O\left(\dfrac x y \log\log x\right)$，空间复杂度为$O(y)$。 1.3 筛法求$\phi(x,a)

咕咕咕...