一份高效的Python埃拉托色尼筛法实现

Question

一份高效的Python埃拉托色尼筛法实现

8

这段非常简短的 #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秒。

有人想分享更有效率的代码以满足上述限制吗？

- Rob

3

像你之前的许多人一样，你试图编写一个筛法，结果却用了试除法。 - user2357112

1

@Roberto，筛法的思想是避免除法，而你的代码并没有这样做。（一个真正的筛法甚至可以没有任何比较就能写出更短的代码。） - MB-F

1

是的，我想指出在这里使用numpy不是一个很好的性能选择。正如kazemakase所指出的，这并不是真正的筛法，因为你要迭代地检查模数。此外，如果您真的想要提高性能，(a!=a[j]) & (a%a[j] == 0)是一种非常浪费的操作。 - juanpa.arrivillaga

2

这里有很多带基准测试的NumPy示例：https://dev59.com/snI95IYBdhLWcg3w-DH0#2068548 - Alex Riley

1

请特别参考这里的答案：https://dev59.com/snI95IYBdhLWcg3w-DH0#3035188 - Alex Riley

显示剩余6条评论

4个回答

3

我决定尝试一下，并创建了一个更优化的 NumPy 版本，由@user2357112 supports 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

在我的笔记本电脑上，我筛选出了前10亿个（1e9）自然数中的质数，用时 0:00:10.378686 秒。JIT编译器提供了至少一个数量级的性能；在撰写本文时，下一个最快的答案需要 0:01:27.059963 分钟。不幸的是，我在这个系统（Mac）上没有Nvidia GPU和Cuda，否则我会使用它们。

- Greenstick

Nvidia Tesla V100 -> 11.798239 秒 - Devashish Das

1

输入的大小是多少？你使用了哪个JIT/参数？ - Greenstick

我使用了你的答案中相同的代码。但是由于某些原因，它没有在nvidia-smi上显示，也没有使用两个GPU。有没有办法让它在多个GPU上运行？ - Devashish Das

1

@DevashishDas 我会查看Numba文档，他们有一个特定的装饰器来利用CUDA。 - Greenstick

1

非常感谢你。使用我的算法需要大约4天才能获得结果。而你的算法只需要3分钟。真的非常感谢你，非常感谢你。 - Srinadh

2

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])

- andreis11

应该注意那里有一个错别字，在推导式中数字变量应该是“limit”。 - Justin Malloy

@JustinMalloy 非常感谢！ - andreis11

1

这里有一个使用NumPy的简单方法来做到这一点。它与 OP 的索引思路相似，而不是在主循环中再次循环，但没有除法检查，只有切片。

它也类似于 user2357112 支持 Monica 的答案，但这个答案只考虑奇数，这使得它更快。

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

1. 创建一个大小为n的布尔数组。 2. 循环遍历小于n的平方根的奇数。
注意:数字可以互换使用作为索引。 3. 将合成数(素数的倍数)标记为false。

最后我们有一个布尔数组，我们可以用它来检查一个小于n的数是否为质数(除了偶数，你可以在不使用数组的情况下检查它们，使用 & 1 或其他方法)。

n = 20000000的平均时间:0.1063秒

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

- Iyad Ahmed

3

以下是有关编程的内容，请将其翻译成中文。请仅返回翻译后的文本：可以包含一些解释，而不是全部代码答案。 - KetZoomer

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user2357112 · Accepted Answer

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

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)

它维护了一个“可能是质数”的标志数组，并直接取消与素数倍数对应的标志，无需测试可被当前处理的素数整除的数字，特别是那些不可被整除的数字。

你正在进行试除法，只是通过检测候选除数是否可被整除来判断。即使是试除法的良好实现也需要执行更多操作和更昂贵的操作，而筛法则不然。您的实现甚至比试除法更耗费工作量，因为它考虑了非质数的候选除数，并且因为它仍在对已知为质数的数字执行可除性测试。