Python 中的高效埃拉托色尼筛法

Question

#Python 中的这段非常简短的代码试图模拟前 N 个自然数的“埃拉托色尼筛法”，并具有 (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 上，它返回第一个中的素数：

• N = 100,000 在 0.03 秒内；
• N = 1,000,000 在 0.63 秒内；
• N = 10,000,000 在 22.2 秒内。

是否有人愿意在上述限制内共享 CPU 时间方面更高效的代码？

Answer 1

Eratosthenes 的一个实际 NumPy 筛子看起来像这样：

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)

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

您正在做的是试用除法，您只需通过并测试数字是否可以被候选除数整除。即使是一个很好的试除法，也需要比筛子做更多的操作，也需要更昂贵的操作。您的实现做的工作比这还多，因为它考虑了非素数候选除数，并且因为它一直在对应该已经知道是素数的数字执行可除性测试。

Answer 2

我决定尝试一下，并创建了一个由 @user2357112 发布的进一步优化的 NumPy 版本，它支持 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 分钟。遗憾的是，我在这个系统 (Mac) 上没有 Nvidia GPU 和 Cuda，否则我会使用它。

Answer 3

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 4

这是使用 NumPy 执行此操作的简单方法。它有点类似于 OP 的索引而不是在主循环内再次循环的想法，但没有除法检查，只有切片。

类似user2357112也支持莫妮卡回答，不过这个只考虑奇数，速度比较快。

这是一个典型的奇数筛：https://stackoverflow.com/a/20248491/8094047。

1. 创建一个大小为 n 的 bool 数组。
2. 循环遍历奇数直到 n 的平方根。
   注意：数字可以互换用作索引
3. 将复合（质数的倍数）数标记为假。

最后我们有一个 bool 数组，我们可以使用它来检查 n 以下的数字是否为素数（偶数除外，您可以在没有数组的情况下检查它们，使用 & 1 左右）

平均。 n = 20000000 的时间：0.1063s

import numpy as np
n = 20000000

isprime = np.ones(n, dtype=np.bool)

# odd only sieve
i = 3
while (i * i < n):
    if isprime[i]:
        isprime[i * i:n:2 * i] = False
    i += 2

# test
print(isprime[97]) # should be True