这是我能想到的最好的算法。
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 time import time
# 175 ms for all the primes up to the value 10**6
def primes_sieve(limit):
a = [True] * limit
a[0] = a[1] = False
#a[2] = True
for n in xrange(4, limit, 2):
a[n] = False
root_limit = int(limit**.5)+1
for i in xrange(3,root_limit):
if a[i]:
for n in xrange(i*i, limit, 2*i):
a[n] = False
return a
LIMIT = 10**6
s=time()
primes = primes_sieve(LIMIT)
print time()-s
下面是我在Python中通常用来生成质数的代码:
$ python -mtimeit -s'import sieve' 'sieve.sieve(1000000)'
10 loops, best of 3: 445 msec per loop
$ cat sieve.py
from math import sqrt
def sieve(size):
prime=[True]*size
rng=xrange
limit=int(sqrt(size))
for i in rng(3,limit+1,+2):
if prime[i]:
prime[i*i::+i]=[False]*len(prime[i*i::+i])
return [2]+[i for i in rng(3,size,+2) if prime[i]]
if __name__=='__main__':
print sieve(100)
它不能与这里发布的更快的解决方案竞争，但至少它是纯python。
谢谢你提出这个问题。我今天真的学到了很多东西。
很抱歉打扰，但erat2()在算法中有一个严重的缺陷。
在搜索下一个合成时，我们只需要测试奇数。
Q p都是奇数;那么q+p是偶数，不需要检验，但q+2*p总是奇数。这消除了while循环条件中的“if even”测试，并节省了大约30%的运行时。
当我们在它:而不是优雅的'D.pop(q,None)'获取和删除方法，使用'if q in D: p=D[q]，del D[q]'，这是两倍的速度!至少在我的机器上(P3-1Ghz)。
所以我建议这个聪明算法的实现:
def erat3( ):
from itertools import islice, count
# q is the running integer that's checked for primeness.
# yield 2 and no other even number thereafter
yield 2
D = {}
# no need to mark D[4] as we will test odd numbers only
for q in islice(count(3),0,None,2):
if q in D: # is composite
p = D[q]
del D[q]
# q is composite. p=D[q] is the first prime that
# divides it. Since we've reached q, we no longer
# need it in the map, but we'll mark the next
# multiple of its witnesses to prepare for larger
# numbers.
x = q + p+p # next odd(!) multiple
while x in D: # skip composites
x += p+p
D[x] = p
else: # is prime
# q is a new prime.
# Yield it and mark its first multiple that isn't
# already marked in previous iterations.
D[q*q] = q
yield q
