Python中Eratosthenes的高效筛选

Question

#Python中的这段非常简短的代码试图模拟前N个自然数的“Eratosthennes筛”，并限制（0）脚本长度；（1） 最小化“if语句”和“for/while循环”；（2） CPU时间方面的效率。

import numpy as np
N = 10**5
a = np.array(range(3,N,2))
for j in range(0, int(round(np.sqrt(N),0))):
    a[(a!=a[j]) & (a%a[j] == 0)] = 0
    a = a[a!=0]
a = [2]+list(a)

在Intel Core I5上，它返回第一个数字中的质数：

< Li > 0.03秒内N = 100，000； < Li > 0.63秒内N = 1，000，000； < li>N = 10，000，000在22.2秒内。

有人想在上述限制内分享更高效的CPU时间代码吗？

Answer 1

10.000.000为1.94秒

def sieve_eratosthene(limit):

    primes = [True] * (limit+1)
    iter = 0

    while iter < limit**0.5 :
        if iter < 2:
            primes[iter]= False

        elif primes[iter]:
            for i in range(iter*2, limit+1, iter):
                primes[i] = False

        iter+=1

    return(x for x in range(limit+1) if primes[x])

Answer 2

我决定玩这个，并创建了一个进一步优化的NumPy版本，由@user2357112发布，支持Monica，它使用Numba JIT来加快速度。

import numba
import numpy
import timeit
import datetime 

@numba.jit(nopython = True, parallel = True, fastmath = True, forceobj = False)
def sieve (n: int) -> numpy.ndarray:
    primes = numpy.full(n, True)
    primes[0], primes[1] = False, False
    for i in numba.prange(2, int(numpy.sqrt(n) + 1)):
        if primes[i]:
            primes[i*i::i] = False
    return numpy.flatnonzero(primes)

if __name__ == "__main__":
    
    timestart = timeit.default_timer()
    print(sieve(1000000000))
    timestop = timeit.default_timer()
    timedelta = (timestop - timestart)
    print(f"Time Elapsed: {datetime.timedelta(seconds = timedelta)}")

else:
    pass

在我的笔记本电脑上，我筛选出0:00:10.378686秒内第一个10亿（1e9）自然数中的素数。JIT在这里提供了至少一个数量级的性能；在撰写本文时，下一个最快的答案需要0:01:27.059963分钟。遗憾的是，我在这个系统（苹果）上没有英伟达图形处理器和库达，否则我会用它。

Answer 3

一个实际的厄拉多塞数字筛看起来像这样:

def sieve(n):
    flags = numpy.ones(n, dtype=bool)
    flags[0] = flags[1] = False
    for i in range(2, n):
        # We could use a lower upper bound for this loop, but I don't want to bother with
        # getting the rounding right on the sqrt handling.
        if flags[i]:
            flags[i*i::i] = False
    return numpy.flatnonzero(flags)

它维护一个“可能素数”标志的数组，并直接取消设置对应于素数倍数的标志，而无需测试可除性，特别是对于当前正在处理的素数不可整除的数字。

我们所做的是试除法，在这里，我们只需要通过并测试数字是否可以被候选除数整除。即使是一个好的试除法实现也需要做更多的操作和更昂贵的操作，而不是筛子。你的实现做的工作甚至更多，因为它考虑了非质数候选除数，并且因为它一直在对它应该知道是质数的数字进行可整除性测试。