这是我能想到的最好的算法。
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
下面是Eratosthenes的一个numpy版本，具有良好的复杂度(低于排序长度为n的数组)和向量化。与@unutbu相比，用46微秒就可以找到100万以下的所有质数。
import numpy as np
def generate_primes(n):
is_prime = np.ones(n+1,dtype=bool)
is_prime[0:2] = False
for i in range(int(n**0.5)+1):
if is_prime[i]:
is_prime[i**2::i]=False
return np.where(is_prime)[0]
计时:
import time
for i in range(2,10):
timer =time.time()
generate_primes(10**i)
print('n = 10^',i,' time =', round(time.time()-timer,6))
>> n = 10^ 2 time = 5.6e-05
>> n = 10^ 3 time = 6.4e-05
>> n = 10^ 4 time = 0.000114
>> n = 10^ 5 time = 0.000593
>> n = 10^ 6 time = 0.00467
>> n = 10^ 7 time = 0.177758
>> n = 10^ 8 time = 1.701312
>> n = 10^ 9 time = 19.322478
这些都是经过编写和测试的。所以没有必要重新发明轮子。
python -m timeit -r10 -s"from sympy import sieve" "primes = list(sieve.primerange(1, 10**6))"
打破了12.2秒的记录!
10 loops, best of 10: 12.2 msec per loop
如果这还不够快，你可以试试PyPy:
pypy -m timeit -r10 -s"from sympy import sieve" "primes = list(sieve.primerange(1, 10**6))"
结果是:
10 loops, best of 10: 2.03 msec per loop
得到247张赞成票的答案列出了15.9毫秒的最佳解决方案。
比较这个! !
下面是Eratosthenes的一个numpy版本，具有良好的复杂度(低于排序长度为n的数组)和向量化。与@unutbu相比，用46微秒就可以找到100万以下的所有质数。
import numpy as np
def generate_primes(n):
is_prime = np.ones(n+1,dtype=bool)
is_prime[0:2] = False
for i in range(int(n**0.5)+1):
if is_prime[i]:
is_prime[i**2::i]=False
return np.where(is_prime)[0]
计时:
import time
for i in range(2,10):
timer =time.time()
generate_primes(10**i)
print('n = 10^',i,' time =', round(time.time()-timer,6))
>> n = 10^ 2 time = 5.6e-05
>> n = 10^ 3 time = 6.4e-05
>> n = 10^ 4 time = 0.000114
>> n = 10^ 5 time = 0.000593
>> n = 10^ 6 time = 0.00467
>> n = 10^ 7 time = 0.177758
>> n = 10^ 8 time = 1.701312
>> n = 10^ 9 time = 19.322478
对于最快的代码，numpy解决方案是最好的。不过，出于纯粹的学术原因，我发布了我的纯python版本，它比上面发布的食谱版本快不到50%。由于我将整个列表放在内存中，所以需要足够的空间来容纳所有内容，但它的可伸缩性似乎相当好。
def daniel_sieve_2(maxNumber):
"""
Given a number, returns all numbers less than or equal to
that number which are prime.
"""
allNumbers = range(3, maxNumber+1, 2)
for mIndex, number in enumerate(xrange(3, maxNumber+1, 2)):
if allNumbers[mIndex] == 0:
continue
# now set all multiples to 0
for index in xrange(mIndex+number, (maxNumber-3)/2+1, number):
allNumbers[index] = 0
return [2] + filter(lambda n: n!=0, allNumbers)
结果是:
>>>mine = timeit.Timer("daniel_sieve_2(1000000)",
... "from sieves import daniel_sieve_2")
>>>prev = timeit.Timer("get_primes_erat(1000000)",
... "from sieves import get_primes_erat")
>>>print "Mine: {0:0.4f} ms".format(min(mine.repeat(3, 1))*1000)
Mine: 428.9446 ms
>>>print "Previous Best {0:0.4f} ms".format(min(prev.repeat(3, 1))*1000)
Previous Best 621.3581 ms
从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
