一秒钟统计出一亿内质数的个数，论Python统计质数的算法优化再优化

关于判断自然数是否质数的算法，网上列举出几种算法，而且速度越来越快，常见的比较快的算法是：“任何一个自然数，总可以表示成以下六种形式之一：6n，6n+1，6n+2，6n+3，6n+4，6n+5（n=0,1,2...）我们可以发现，除了2和3，只有形如6n+1和6n+5的数有可能是质数。且形如6n+1和6n+5的数如果不是质数，它们的因数也会含有形如6n+1或者6n+5的数”，因此用Python可以得到如下算法代码：

方法1：

# 统计一亿以内的质数
import time
from numba import jit

# 方法1
# 引入numba库的jit装饰器进行加速运算
@jit
def test1(end):
    count = 0
    start = 0
    def step_two(x):
        for i in range(5, int(x ** 0.5) + 1, 6):
            if (x % i == 0) or (x % (i + 2) == 0):
                return False
        return True
        
    if end >3:
        count = 2
        start = 4
    for x in range(start,end):
        if (x % 6 != 1) and (x % 6 != 5):
            continue
        else:
            if step_two(x):
                count +=1
    return count

if __name__ == '__main__':
    end:int = 100000000
    now = time.time()
    count = test1(end)
    print(f'方法1：\n{end:,}以内的质数有{count}个\n耗时：{time.time()-now:.3f}秒。')

为加快运算，该代码已引入numba库的jit装饰器，统计一亿内的质数个数，在我的台式机上运行耗时84秒。如果没有这个装饰器，耗时会延长几倍。

方法2

我在外网找到更优化的算法，被称为Sieve算法，先创建一个长度为n、只有True值的长列表，然后把所有公倍数的元素值设为False并排除，最后剩下得到的True元素皆为质数。

参考来源：

Analysis of Different Methods to find Prime Number in Python

# 统计一亿以内的质数

import time
from numba import jit

# 引入numba库的jit装饰器进行加速运算
@jit
def SieveOfEratosthenes(n):
    prime = [True for i in range(n+1)]
    p = 2
    pp = p * p
    while(pp <= n):
        if (prime[p] == True):
            for i in range(pp, n + 1, p):
                prime[i] = False
        p += 1
        pp = p * p
    count = 0
    for p in range(2, n):
        if prime[p]:
            count += 1
    return count

if __name__ == '__main__':
    end:int = 100000000
    # 方法2
    now = time.time()
    count = SieveOfEratosthenes(end)
    print(f'方法2：\n{end:,}以内的质数有{count}个\n耗时：{time.time()-now:.3f}秒。')

Sieve算法只有正整数在运算，没有浮点运算，在jit装饰器的加速作用下快得惊人，统计一亿内的质数个数仅需2.084秒！速度比方法1提升35倍！

这算法的速度到天花板了吗？还没有！

方法3

引入Numpy库可以再加速运算。根据方法2，改用Numpy库的语句优化一下代码。

（1）

'''
prime = [True for i in range(n+1)]
这一句是创建 n 个True值的数组，使用numpy语句可改为：
'''
prime = np.full(n, True)

（2）

'''
for i in range(pp, n + 1, p):
    prime[i] = False

遍历循环可以浓缩成一行代码：
'''
prime[pp::p] = False

（3）

'''
count = 0
for p in range(2, n):
    if prime[p]:
        count += 1
return count

统计prime数组里True值的个数。
使用numpy语句浓缩成一行代码：
'''
return np.count_nonzero(prime)-2

这里return的后面为什么要减2？是因为在统计质数的算法中，0和1从严格意义上来说不算是质数，prime数组的第0、第1个元素就是对应0和1，为了与方法1、方法2的计算结果保持一致，故从结果中减去2。

方法3的完整代码如下：

from time import time
import numpy as np

def test3(n):
    prime = np.full(n, True)
    p = 2
    pp = p * p
    while(pp <= n):
        if (prime[p] == True):
            prime[pp::p] = False
        p += 1
        pp = p * p
    count=np.count_nonzero(prime) - 2
    return count

if __name__=='__main__':
    end:int = 100000000
    now = time()
    count = test3(end)
    print(f'方法3：\n{end:,}以内的质数有{count}个\n耗时：{time()-now:.3f}秒。')

在我主频为3GHz的台式电脑上运行结果仅为0.963秒！运算速度比方法2提升一倍，比方法1提升87倍！

 

注意到什么没有？

方法3的使用了numpy情况下，不需要加上@jit的装饰器。

因为numpy库的运算速度本来就很快。如果加上了@jit，反而拖累了运算速度。