这是我能想到的最好的算法。
def get_primes(n):
numbers = set(range(n, 1, -1))
primes = []
while numbers:
p = numbers.pop()
primes.append(p)
numbers.difference_update(set(range(p*2, n+1, p)))
return primes
>>> timeit.Timer(stmt='get_primes.get_primes(1000000)', setup='import get_primes').timeit(1)
1.1499958793645562
还能做得更快吗?
这段代码有一个缺陷:由于numbers是一个无序集,不能保证numbers.pop()将从集合中移除最低的数字。尽管如此,它还是适用于(至少对我来说)一些输入数字:
>>> sum(get_primes(2000000))
142913828922L
#That's the correct sum of all numbers below 2 million
>>> 529 in get_primes(1000)
False
>>> 529 in get_primes(530)
True
这是问题解的一种变化应该比问题本身更快。它使用埃拉托色尼的静态筛,没有其他优化。
from typing import List
def list_primes(limit: int) -> List[int]:
primes = set(range(2, limit + 1))
for i in range(2, limit + 1):
if i in primes:
primes.difference_update(set(list(range(i, limit + 1, i))[1:]))
return sorted(primes)
>>> list_primes(100)
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
从2021年的答案开始,我还没有发现二进制数组方法对10亿以下的质数有利。
但我可以用几个技巧将质数从2加速到接近x2:
使用numexpr库将numpy表达式转换为分配较少的紧循环
取代np。有更快的选择
以某种方式操作筛选的前9个元素,因此不需要改变数组的形状
总之,在我的机器上,质数< 10亿的时间从25秒变成了14.5秒
import numexpr as ne
import numpy as np
def primesfrom2to_numexpr(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Input n>=24, Returns a array of primes, 2 <= p < n + a few over"""
sieve = np.zeros((n // 3 + (n % 6 == 2))//4+1, dtype=np.int32)
ne.evaluate('sieve + 0x01010101', out=sieve)
sieve = sieve.view('int8')
#sieve = np.ones(n // 3 + (n % 6 == 2), dtype=np.bool_)
sieve[0] = 0
for i in np.arange(int(n ** 0.5) // 3 + 1):
if sieve[i]:
k = 3 * i + 1 | 1
sieve[((k * k) // 3)::2 * k] = 0
sieve[(k * k + 4 * k - 2 * k * (i & 1)) // 3::2 * k] = 0
sieve[[0,8]] = 1
result = np.flatnonzero(sieve)
ne.evaluate('result * 3 + 1 + result%2', out=result)
result[:9] = [2,3,5,7,11,13,17,19,23]
return result
使用Numpy实现的半筛子略有不同:
http://rebrained.com/?p=458
import math
import numpy
def prime6(upto):
primes=numpy.arange(3,upto+1,2)
isprime=numpy.ones((upto-1)/2,dtype=bool)
for factor in primes[:int(math.sqrt(upto))]:
if isprime[(factor-2)/2]: isprime[(factor*3-2)/2:(upto-1)/2:factor]=0
return numpy.insert(primes[isprime],0,2)
有人能把这个和其他时间比较一下吗?在我的机器上,它似乎与其他Numpy半筛相当。