在做欧拉计划的习题时，关于素数的问题比较常见。现在的人类数学暂时还没有素数通项公式，只有一个定义，但这不能阻止人类寻找各种关于素数的问题的解决办法。通过素数的定义，直观地使用试除法判断一个数是不是素数，其效率极其低下，时间复杂度为\(O(n^2)\)。但是数学家们在寻找素数的出现规律时，通过统计，发现小于自然数n的范围内，素数出现的概率是有一定的规律的（虽然还是没摸清楚），这就发展出素数定理。
这里不是讨论素数定理。数学家们在求小于n的素数有多少个、它们分别都是什么的时候，发明了筛选法，以便快速求出素数表。这里列举出现时为止常见的素数筛选法。
埃拉托斯特尼筛选法
埃拉托斯特尼（Eratosthenes, 276BC ~ 194BC），古希腊数学家，就是那位人类有历史记载以来第一位估算出地球直径、设立经纬度系统的数学家。他通过列出所有小于n的数，从小至大遍历这n-1个数，然后利用合数的定义，将后面那些明显是较小的数的倍数标记为非素数，如此类推，没被标记为非素数的数，便是所有小于n的素数，整个过程正如标题图片那样，像筛子那样一次次将不要的东西筛走，这个算法堪称筛选式算法的鼻祖。
该算法时间复杂度为\(O(n \log \log n)\)：
// Eratosthenes Sieve
const start = Date.now()
function* primesBelow(n) {
const sieve = new BitSet // 1：是素数，0：非素数
sieve.flip(0, n) // 全部设置为1
.set(0, 0) // 忽略0、1
.set(1, 0)
for (let i = 2; i < n; ++i) {
if (sieve.get(i) == 0)
continue
for (let j = 2 * i; j < n; j += i)
sieve.set(j, 0)
}
for (let i = 2; i < n; ++i)
if (sieve.get(i) == 1)
yield i
}
var sum = 0
for (const prime of primesBelow(5000000))
sum += prime
alert(sum + ', ' + (Date.now() - start) + 'ms')
欧拉筛选法
欧拉（Euler, 1707 ~ 1783），旧瑞士联邦巴塞尔人，是理科生无人不晓的近代数学家。他当时在证明黎曼ζ函数的乘积公式时，发现埃拉托斯特尼筛选法不够快，因为它会重复对两个素数的乘积的倍数进行判断，例如筛走了2的倍数后，在筛选3的倍数时，会再一次碰见6、12、24这些数。于是他便通过记录前面已经确定是素数的数字，消除重复筛走倍数的过程。用计算机算法的术语讲，这个算法是利用空间换时间。
该算法时间复杂度为\(O(n)\)：
// Euler Sieve
const start = Date.now()
const primeCountUpBound = (n) => Math.ceil(n < 55 ? 1.25506*n/Math.log(n) : n/(Math.log(n)-4))
function* primesBelow(n) {
const sieve = new BitSet // 0：是素数，1：非素数
const primes = new Array(primeCountUpBound(n)) // 利用素数定理一次性分配空间，避免数组重复申请空间
var count = 0
for (let i = 2; i < n; ++i) {
if (sieve.get(i) == 0) {
primes[count++] = i // 发现一个素数便记录一个
yield i
}
for (let j = 0; j < count; ++j) {
const prime = primes[j]
if (i * prime >= n)
break
sieve.set(i * prime, 1)
if ((i % prime) == 0)
break // 这一跳是控制不要往太后筛走有可能重复的数字的关键
}
}
}
var sum = 0
for (const prime of primesBelow(5000000))
sum += prime
alert(sum + ', ' + (Date.now() - start) + 'ms')
欧拉筛选法相比埃拉托斯特尼筛选法，最大的差异，是埃拉托斯特尼筛选法是在遍历的过程中，用自然数与当前的数字逐个相乘筛走倍数，而欧拉的方法是使当前数字与小于等于它的素数逐个相乘筛走倍数，所以要额外的空间记录素数。欧拉筛选法这种筛走倍数的策略，不像埃拉托斯特尼筛选法每次都从当前数字到n都筛走一遍，而是只筛走当前数字后面一小段范围内的倍数，随着记录下来的素数越多，它一次性筛走的倍数就越多（但受n的限制后面就变得越来越少），而且还有一个规则，如果当前的数字是一个合数，那么它的素因数肯定能够在已找出来的素数中找到，那么只标记筛走该数字与2至第一个素因数相乘的倍数，即人为限制形如\(m p^2\)的合数只被标记一次。
以少于21的情况为例，加粗、加下划线的为当前选中的数字，加删除线的为被筛走的数字，红色代表当前数字与素数相乘的合数：
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 primes: [2]
2 345 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 primes: [2, 3]
2 34567 8910 4自身能被已找出来的素数2整除，即它自身是个合数，只标记一个倍数，然后跳过 11 12 13 14 15 16 17 18 19 20 primes: [2, 3]
2 345678910 11 12 13 14 15 16 17 18 19 20 primes: [2, 3, 5]
2 3456789106自身能被已找出来的素数2整除，即它自身是个合数，只标记一个倍数，然后跳过 11 12 13 141516 17 18 19 20 primes: [2, 3, 5]
2 345678910111213 141516 17 18 19 20 primes: [2, 3, 5, 7]
2 345678910111213141516 17 18 19 20 primes: [2, 3, 5, 7]
2 34567891011121314151617 18 19 20 primes: [2, 3, 5, 7]
2 345678910111213141516171819 20 primes: [2, 3, 5, 7]
……
2 34567891011121314151617181920primes: [2, 3, 5, 7, 11, 13, 17, 19]
Sundaram筛选法
S.P. Sundaram，近代印度人，具体信息不明，人们只是传说他在1934年发明了这个算法，当时他还是一名学生。该算法是埃拉托斯特尼筛选法的变种，它从奇数的方向进行突破。
观察到所有素数除了2其他都是奇数，那么，先将将奇数表达为：\(2m + 1, m \in \mathbb{N}\)。在这些奇数中，将那些合数排除掉，那么剩下的便都是素数。
假设\(i,j \in \mathbb{N}, 1 < i \le j, m = i + j + 2ij\)，那么形如\(2(i + j + 2ij) + 1\)的奇数一定是合数，因为：
- \(2(i + j + 2ij) + 1\)；
- \(2i + 2j + 4ij + 1\)；
- \((2i + 1)(2j + 1)\)；
- 两个奇数的乘积肯定是一个合数。
所以说，只要从1到m找出能满足\(m = i + j + 2ij\)的所有\((i, j)\)组合，就相当于将范围小于\(n = 2m + 2\)的所有奇合数都筛选出来，剩下的使用奇数表达式就能算出素数。
相当简单易明但又极其天才的算法，虽然表面上和欧拉筛选法的时间复杂度都是\(O(n)\)，但是它有效地将整个筛选器的大小缩小了一半，还不需要占用额外的空间记录已找出来的素数：
// Sundaram Sieve
const start = Date.now()
function* primesBelow(n) {
const m = Math.floor((n - 2) / 2) // 从n = 2m + 2反推
const sieve = new BitSet // 0：保留，1：筛走
for (let i = 1; i < m; ++i) {
let j = i, p = i + j + 2*i*j
for (; p <= m; ++j, p = i + j + 2*i*j)
sieve.set(p, 1) // 没超范围的只要满足上述公式的都筛走
}
yield 2
for (let i = 1; i <= m; ++i)
if (sieve.get(i) == 0)
yield (2*i + 1) // 重新组合计算出素数
}
var sum = 0
for (const prime of primesBelow(5000000))
sum += prime
alert(sum + ', ' + (Date.now() - start) + 'ms')
印度就是这么神奇，能想出这个算法的人竟然是个学生，而且到后来也不知道他遭遇了什么，一直籍籍无名。拉马努金是比他早一两辈的人（一个死了另一个才出生），要不是机缘巧合去了英国，他搞不好也一辈子在印度教庙宇里写粉笔。
Atkin筛选法
A. O. L. Atkin（1925 ~ 2008），英国人，在2003年与Daniel J. Bernstein一起发明了这个算法。这个算法比较复杂难懂，使用了好几种条件进行取反操作，还做了形如\(kn^2\)的合数的快速筛选。心情差一点、天气热一点都没耐心看完这个算法到底是如何工作、证明的。?
该算法时间复杂度为\(O(n)\)，实际执行起来效率和欧拉筛选法差不多，但是却不需要占用额外的空间记录已找出来的素数：
// Atkin Sieve
const start = Date.now()
function* primesBelow(limit) {
const sqrtLimit = Math.floor(Math.sqrt(limit))
const sieve = new BitSet // 0：非素数，1：素数
sieve.set(0, 0)
.set(1, 0)
.set(2, 1)
.set(3, 1)
for (let x = 1; x <= sqrtLimit; ++x) {
const xx = x*x
for (let y = 1; y <= sqrtLimit; ++y) {
const yy = y*y
let n = 4*xx + yy
if (n <= limit && ((n % 12) == 1 || (n % 12) == 5))
sieve.flip(n)
n = 3*xx + yy
if (n <= limit && (n % 12) == 7)
sieve.flip(n)
n = 3*xx - yy
if (x > y && n <= limit && (n % 12) == 11)
sieve.flip(n)
}
}
for (let n = 5; n <= sqrtLimit; ++n)
if (sieve.get(n) == 1) {
const nn = n*n
for (let p = nn; p <= limit; p += nn)
sieve.set(p, 0)
}
for (let i = 0; i <= limit; ++i)
if (sieve.get(i) == 1)
yield i
}
var sum = 0
for (const prime of primesBelow(5000000))
sum += prime
alert(sum + ', ' + (Date.now() - start) + 'ms')
算法解释：
……未完待续（To Be Continued）……
效率比较
耗时：Sundaram < 欧拉 = Atkin < 埃拉托斯特尼。
空间占用：Sundaram < 埃拉托斯特尼 = Atkin < 欧拉。
即便连上算法理解的难易程度，Sundaram筛选器都是相当出色的算法。
打赏作者
由 bruce 于 2020-06-8 19:22 更新
