技术

如何快速查找质数

质数（素数），指在大于 1 的自然数中，除了 1 和该数自身外，无法被其他自然数整除的数（如 2，3，5，7，11 ……）。那么，如何能简单快速地找出质数呢？

埃拉托斯特尼筛法

埃拉托斯特尼（Eratosthenes）是一位古希腊数学家、地理学家、历史学家、诗人、天文学家，也是阿基米德的好基友。除了发明质数的筛选方法，他还测量地球周长、日地间距、地月间距，编排星图，绘制地图，把名字写在月球上（埃拉托斯特尼陨石坑）。

埃拉托斯特尼筛法（sieve of Eratosthenes ）用来找出一定范围（n）内的所有质数。其方法是从 2 开始，在$\sqrt{n}$ 以内，将每个质数的倍数剔除掉，剩下的就是所求范围的质数。例如找 100 以内的质数，先把 2 的倍数筛掉（保留 2），再把 3 的倍数筛掉（保留 3），如此重复下去，直到 7 的倍数被筛掉（因为下一个质数 11 已经大于$\sqrt{100}$ )，剩下的就是 100 以内的质数。


一个亿的小目标

王先生说： “最好先定一个能达到的小目标，比方说我先挣它一个亿”。挣一个亿，我太难了，不如先 “算一个亿”吧。下面我们就用埃拉托斯特尼筛法来找 1 亿以内的质数，为了方便，结果只输出该范围内质数的个数。

Python

  • 方式一：列表

先创建一个包含 N+1 个 True 的列表，然后把非质数设为 False，最后剩下的 True 就是质数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
N = 100000000
primes = [ True for i in range(N) ]
primes[:1] = [ False, False ]  # 去掉0和1

k = 2 
while k * k <= N:
    if primes[k]:
        j = k * k
        while j <= N:
            primes[j] = False
            j += k
    k += 1

print("{} primes in {}".format(primes.count(True), N))

结果：

1
5761455 primes in 100000000

耗时 67 秒。

  • 方式二：集合

我们还可以用集合的方式来筛选，先创建一个从 3 到 N 的奇数集合（加上 2），然后减去每个质数的倍数的集合，最后得到只有质数的集合（代码缩减了一半）。

1
2
3
4
5
6
7
8
N = 100000000
primes = {i for i in range(3, N + 1, 2)} | {2}

for j in range(3, int(N ** 0.5) + 1, 2):
    if j in primes:
        primes -= {i for i in range(j * j, N + 1, j)}

print("{} primes in {}:".format(len(primes), N))

耗时 59 秒，有进步。（ 仿佛听到有人喊： “偷步啊！你从奇数开始筛。“ 答：“是，又怎样😆！” ）

  • 方式三：函数

看到 Wiki 上 埃拉托斯特尼筛法 的 python 代码，写得真简练，拿来试试。

1
2
3
4
5
6
7
8
9
10
11
def eratosthenes(n):
     IsPrime = [True] * (n + 1)
     for i in range(2, int(n ** 0.5) + 1): 
         if IsPrime[i]:
             for j in range(i * i, n + 1, i): 
                 IsPrime[j] = False
     return {x for x in range(2, n + 1) if IsPrime[x]}

N = 100000000
primes = eratosthenes(N)
print("{} primes in {}:".format(len(primes), N))

耗时一下子降到 26 秒，凭什么啊，如果把函数拆了，会增加 10 秒，暂时没搞懂。

Go

现在算法已经不是问题，时间才是问题，我要尽快算到一个亿。于是我写了一段憋足的 Go 代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
package main
import "fmt"

func main() {
    N := 100000000
    primes := make([]bool, N+1)
    for i := 2; i <= N; i += 1 { 
        primes[i] = true
    }   

    for k := 2; k*k <= N; k += 1 { 
        if primes[k] {
           for j := k*k; j <= N; j += k { 
               primes[j] = false
           }   
       }   
   }   

   total := 0
   for i := 0; i <= N; i++ {
       if primes[i]{
           total += 1
       }   
   }   
   fmt.Println(total, "primes in", N)

}

go run 一下，1 个亿耗时仅 1.5 秒！这速度堪比双十一的第 1 秒。

在一个亿面前，Python 被 Go 打得满地找牙，解释型就是打不过编译型吗？

这时，Python 背后传来一个坚实的声音：“别怕，有我在！”

Python + C++

这个背后的声音是 primesieve-python，一个调用 primesieve C++ library 的 Python 模块。

安装： pip3 install primesieve

少啰嗦，给我上！

1
2
3
4
import primesieve
N = 100000000
prime = primesieve.primes(N)
print(len(prime), "primes in", N)

1 个亿耗时 0.5 秒 ！10 个亿 4.5 秒 ！

如果只需要质数的个数，直接调用 primesieve.count_primes(N)，1 亿只需 0.06 秒，1000 亿只需 30 秒。

Python 成功逆袭，Go 表示请外援不算好汉～


总结

最后，我们统计一下以上各种求质数方式的效率（耗时）：

* 系统环境：MacOS 10.14.4 / 2.6 GHz Intel 2-Core i5 / 8GB DDR3 / SSD

N以内 质数的个数 Python / 列表 Python / 集合 Python / 函数 Go Python / primesieve
10万 9592 0.09s 0.07s 0.10s 0.006s 0.03s
100万 78498 0.6s 0.3s 0.25s 0.01s 0.04s
1000万 664579 6.1s 3.7s 2.6s 0.10s 0.11s
1亿 5761455 1m7s 54s 26s 1.6s 0.54s
10亿 50847534 算了吧 压力太大 太难了 15s 4.5s
20亿 98222287 32s 10s
30亿 144449537 53s 24s
40亿 189961812 1m17s 44s
50亿 234954223 1m49s 1m1s
60亿 279545368 2m31s 1m25s
100亿 455052511 >100m 以下用 count_primes(n)
2s
200亿 882206716 参考楼上 5s
300亿 1300005926 8s
400亿 1711955433 11s
500亿 2119654578 14s
1000亿 4118054813 30s

（完）