这是我能想到的最好的算法。
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
这里有一个来自Python Cookbook的非常简洁的示例——该URL的最快版本是:
import itertools
def erat2( ):
D = { }
yield 2
for q in itertools.islice(itertools.count(3), 0, None, 2):
p = D.pop(q, None)
if p is None:
D[q*q] = q
yield q
else:
x = p + q
while x in D or not (x&1):
x += p
D[x] = p
这就给出了
def get_primes_erat(n):
return list(itertools.takewhile(lambda p: p<n, erat2()))
在shell提示符(正如我喜欢做的那样)中测量这段代码在pri.py中，我观察到:
$ python2.5 -mtimeit -s'import pri' 'pri.get_primes(1000000)'
10 loops, best of 3: 1.69 sec per loop
$ python2.5 -mtimeit -s'import pri' 'pri.get_primes_erat(1000000)'
10 loops, best of 3: 673 msec per loop
所以看起来食谱解决方案的速度是原来的两倍多。
很抱歉打扰，但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
随着时间的推移，我收集了几个质数筛子。我电脑上最快的是这样的:
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
我可能迟到了，但必须为此添加自己的代码。它使用大约n/2的空间，因为我们不需要存储偶数，我还使用bitarray python模块，进一步大幅减少内存消耗，并允许计算所有高达1,000,000,000的质数
from bitarray import bitarray
def primes_to(n):
size = n//2
sieve = bitarray(size)
sieve.setall(1)
limit = int(n**0.5)
for i in range(1,limit):
if sieve[i]:
val = 2*i+1
sieve[(i+i*val)::val] = 0
return [2] + [2*i+1 for i, v in enumerate(sieve) if v and i > 0]
python -m timeit -n10 -s "import euler" "euler.primes_to(1000000000)"
10 loops, best of 3: 46.5 sec per loop
这是在64bit 2.4GHZ MAC OSX 10.8.3上运行的
使用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半筛相当。
