这是我能想到的最好的算法。
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 3
def rwh_primes2(n):
correction = (n%6>1)
n = {0:n,1:n-1,2:n+4,3:n+3,4:n+2,5:n+1}[n%6]
sieve = [True] * (n//3)
sieve[0] = False
for i in range(int(n**0.5)//3+1):
if sieve[i]:
k=3*i+1|1
sieve[ ((k*k)//3) ::2*k]=[False]*((n//6-(k*k)//6-1)//k+1)
sieve[(k*k+4*k-2*k*(i&1))//3::2*k]=[False]*((n//6-(k*k+4*k-2*k*(i&1))//6-1)//k+1)
return [2,3] + [3*i+1|1 for i in range(1,n//3-correction) if sieve[i]]
很抱歉打扰,但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
假设N < 9,080,191, Miller-Rabin's Primality检验的确定性实现
import sys
def miller_rabin_pass(a, n):
d = n - 1
s = 0
while d % 2 == 0:
d >>= 1
s += 1
a_to_power = pow(a, d, n)
if a_to_power == 1:
return True
for i in range(s-1):
if a_to_power == n - 1:
return True
a_to_power = (a_to_power * a_to_power) % n
return a_to_power == n - 1
def miller_rabin(n):
if n <= 2:
return n == 2
if n < 2_047:
return miller_rabin_pass(2, n)
return all(miller_rabin_pass(a, n) for a in (31, 73))
n = int(sys.argv[1])
primes = [2]
for p in range(3,n,2):
if miller_rabin(p):
primes.append(p)
print len(primes)
根据维基百科(http://en.wikipedia.org/wiki/Miller -Rabin_primality_test)上的文章,对于a = 37和73,测试N < 9,080,191足以判断N是否为合数。
我从原始米勒-拉宾测试的概率实现中改编了源代码:https://www.literateprograms.org/miller-rabin_primality_test__python_.html
随着时间的推移,我收集了几个质数筛子。我电脑上最快的是这样的:
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 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
所以看起来食谱解决方案的速度是原来的两倍多。
在Pure Python中最快的质数筛分:
from itertools import compress
def half_sieve(n):
"""
Returns a list of prime numbers less than `n`.
"""
if n <= 2:
return []
sieve = bytearray([True]) * (n // 2)
for i in range(3, int(n ** 0.5) + 1, 2):
if sieve[i // 2]:
sieve[i * i // 2::i] = bytearray((n - i * i - 1) // (2 * i) + 1)
primes = list(compress(range(1, n, 2), sieve))
primes[0] = 2
return primes
我优化了埃拉托色尼筛子的速度和内存。
基准
from time import clock
import platform
def benchmark(iterations, limit):
start = clock()
for x in range(iterations):
half_sieve(limit)
end = clock() - start
print(f'{end/iterations:.4f} seconds for primes < {limit}')
if __name__ == '__main__':
print(platform.python_version())
print(platform.platform())
print(platform.processor())
it = 10
for pw in range(4, 9):
benchmark(it, 10**pw)
输出
>>> 3.6.7
>>> Windows-10-10.0.17763-SP0
>>> Intel64 Family 6 Model 78 Stepping 3, GenuineIntel
>>> 0.0003 seconds for primes < 10000
>>> 0.0021 seconds for primes < 100000
>>> 0.0204 seconds for primes < 1000000
>>> 0.2389 seconds for primes < 10000000
>>> 2.6702 seconds for primes < 100000000