在做欧拉计划的习题时,关于素数的问题比较常见。现在的人类数学暂时还没有素数通项公式,只有一个定义,但这不能阻止人类寻找各种关于素数的问题的解决办法。通过素数的定义,直观地使用试除法判断一个数是不是素数,其效率极其低下,时间复杂度为\(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 更新