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
#1 楼
警告：由于硬件或Python版本的差异，timeit的结果可能会有所不同。 下面是一个脚本，它比较了许多实现：
rwh_primes，
rwh_primes1，
rwh_primes2，
sieveOfAtkin，
sieveOfEratosthenes，
sundaram3，
sieve_wheel_30，
ambi_sieve（需要numpy）
primesfrom3to（需要numpy）
primesfrom2to（需要numpy）
非常感谢stephan带来了sieve_wheel_30。
贷给罗伯特·威廉·汉克斯（Robert William Hanks）的primesfrom2to，primesfrom3to，rwh_primes，rwh_primes1和rwh_primes2。
rwh_primes1是最快的测试方法。
+---------------------+-------+
| Method | ms |
+---------------------+-------+
| rwh_primes1 | 43.0 |
| sieveOfAtkin | 46.4 |
| rwh_primes | 57.4 |
| sieve_wheel_30 | 63.0 |
| rwh_primes2 | 67.8 |
| sieveOfEratosthenes | 147.0 |
| ambi_sieve_plain | 152.0 |
| sundaram3 | 194.0 |
+---------------------+-------+
在没有psyco的普通Python方法中，n = 1000000，
rwh_primes2是最快的。
+---------------------+-------+
| Method | ms |
+---------------------+-------+
| rwh_primes2 | 68.1 |
| rwh_primes1 | 93.7 |
| rwh_primes | 94.6 |
| sieve_wheel_30 | 97.4 |
| sieveOfEratosthenes | 178.0 |
| ambi_sieve_plain | 286.0 |
| sieveOfAtkin | 314.0 |
| sundaram3 | 416.0 |
+---------------------+-------+
测试的所有方法中，允许小数y，对于n = 1000000，
primesfrom2to是测试最快的。
+---------------------+-------+
| Method | ms |
+---------------------+-------+
| primesfrom2to | 15.9 |
| primesfrom3to | 18.4 |
| ambi_sieve | 29.3 |
+---------------------+-------+
使用以下命令测量时间：
python -mtimeit -s"import primes" "primes.{method}(1000000)"
将
{method}替换为每个方法名称。#!/usr/bin/env python
import psyco; psyco.full()
from math import sqrt, ceil
import numpy as np
def rwh_primes(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Returns a list of primes < n """
sieve = [True] * n
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i]:
sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
return [2] + [i for i in xrange(3,n,2) if sieve[i]]
def rwh_primes1(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Returns a list of primes < n """
sieve = [True] * (n/2)
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i/2]:
sieve[i*i/2::i] = [False] * ((n-i*i-1)/(2*i)+1)
return [2] + [2*i+1 for i in xrange(1,n/2) if sieve[i]]
def rwh_primes2(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Input n>=6, Returns a list of primes, 2 <= p < 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 xrange(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 xrange(1,n/3-correction) if sieve[i]]
def sieve_wheel_30(N):
# http://zerovolt.com/?p=88
''' Returns a list of primes <= N using wheel criterion 2*3*5 = 30
Copyright 2009 by zerovolt.com
This code is free for non-commercial purposes, in which case you can just leave this comment as a credit for my work.
If you need this code for commercial purposes, please contact me by sending an email to: info [at] zerovolt [dot] com.'''
__smallp = ( 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, 101, 103, 107, 109, 113, 127, 131, 137, 139,
149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311,
313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401,
409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491,
499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683,
691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887,
907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997)
wheel = (2, 3, 5)
const = 30
if N < 2:
return []
if N <= const:
pos = 0
while __smallp[pos] <= N:
pos += 1
return list(__smallp[:pos])
# make the offsets list
offsets = (7, 11, 13, 17, 19, 23, 29, 1)
# prepare the list
p = [2, 3, 5]
dim = 2 + N // const
tk1 = [True] * dim
tk7 = [True] * dim
tk11 = [True] * dim
tk13 = [True] * dim
tk17 = [True] * dim
tk19 = [True] * dim
tk23 = [True] * dim
tk29 = [True] * dim
tk1[0] = False
# help dictionary d
# d[a , b] = c ==> if I want to find the smallest useful multiple of (30*pos)+a
# on tkc, then I need the index given by the product of [(30*pos)+a][(30*pos)+b]
# in general. If b < a, I need [(30*pos)+a][(30*(pos+1))+b]
d = {}
for x in offsets:
for y in offsets:
res = (x*y) % const
if res in offsets:
d[(x, res)] = y
# another help dictionary: gives tkx calling tmptk[x]
tmptk = {1:tk1, 7:tk7, 11:tk11, 13:tk13, 17:tk17, 19:tk19, 23:tk23, 29:tk29}
pos, prime, lastadded, stop = 0, 0, 0, int(ceil(sqrt(N)))
# inner functions definition
def del_mult(tk, start, step):
for k in xrange(start, len(tk), step):
tk[k] = False
# end of inner functions definition
cpos = const * pos
while prime < stop:
# 30k + 7
if tk7[pos]:
prime = cpos + 7
p.append(prime)
lastadded = 7
for off in offsets:
tmp = d[(7, off)]
start = (pos + prime) if off == 7 else (prime * (const * (pos + 1 if tmp < 7 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 11
if tk11[pos]:
prime = cpos + 11
p.append(prime)
lastadded = 11
for off in offsets:
tmp = d[(11, off)]
start = (pos + prime) if off == 11 else (prime * (const * (pos + 1 if tmp < 11 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 13
if tk13[pos]:
prime = cpos + 13
p.append(prime)
lastadded = 13
for off in offsets:
tmp = d[(13, off)]
start = (pos + prime) if off == 13 else (prime * (const * (pos + 1 if tmp < 13 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 17
if tk17[pos]:
prime = cpos + 17
p.append(prime)
lastadded = 17
for off in offsets:
tmp = d[(17, off)]
start = (pos + prime) if off == 17 else (prime * (const * (pos + 1 if tmp < 17 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 19
if tk19[pos]:
prime = cpos + 19
p.append(prime)
lastadded = 19
for off in offsets:
tmp = d[(19, off)]
start = (pos + prime) if off == 19 else (prime * (const * (pos + 1 if tmp < 19 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 23
if tk23[pos]:
prime = cpos + 23
p.append(prime)
lastadded = 23
for off in offsets:
tmp = d[(23, off)]
start = (pos + prime) if off == 23 else (prime * (const * (pos + 1 if tmp < 23 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 29
if tk29[pos]:
prime = cpos + 29
p.append(prime)
lastadded = 29
for off in offsets:
tmp = d[(29, off)]
start = (pos + prime) if off == 29 else (prime * (const * (pos + 1 if tmp < 29 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# now we go back to top tk1, so we need to increase pos by 1
pos += 1
cpos = const * pos
# 30k + 1
if tk1[pos]:
prime = cpos + 1
p.append(prime)
lastadded = 1
for off in offsets:
tmp = d[(1, off)]
start = (pos + prime) if off == 1 else (prime * (const * pos + tmp) )//const
del_mult(tmptk[off], start, prime)
# time to add remaining primes
# if lastadded == 1, remove last element and start adding them from tk1
# this way we don't need an "if" within the last while
if lastadded == 1:
p.pop()
# now complete for every other possible prime
while pos < len(tk1):
cpos = const * pos
if tk1[pos]: p.append(cpos + 1)
if tk7[pos]: p.append(cpos + 7)
if tk11[pos]: p.append(cpos + 11)
if tk13[pos]: p.append(cpos + 13)
if tk17[pos]: p.append(cpos + 17)
if tk19[pos]: p.append(cpos + 19)
if tk23[pos]: p.append(cpos + 23)
if tk29[pos]: p.append(cpos + 29)
pos += 1
# remove exceeding if present
pos = len(p) - 1
while p[pos] > N:
pos -= 1
if pos < len(p) - 1:
del p[pos+1:]
# return p list
return p
def sieveOfEratosthenes(n):
"""sieveOfEratosthenes(n): return the list of the primes < n."""
# Code from: <dickinsm@gmail.com>, Nov 30 2006
# http://groups.google.com/group/comp.lang.python/msg/f1f10ced88c68c2d
if n <= 2:
return []
sieve = range(3, n, 2)
top = len(sieve)
for si in sieve:
if si:
bottom = (si*si - 3) // 2
if bottom >= top:
break
sieve[bottom::si] = [0] * -((bottom - top) // si)
return [2] + [el for el in sieve if el]
def sieveOfAtkin(end):
"""sieveOfAtkin(end): return a list of all the prime numbers <end
using the Sieve of Atkin."""
# Code by Steve Krenzel, <Sgk284@gmail.com>, improved
# Code: https://web.archive.org/web/20080324064651/http://krenzel.info/?p=83
# Info: http://en.wikipedia.org/wiki/Sieve_of_Atkin
assert end > 0
lng = ((end-1) // 2)
sieve = [False] * (lng + 1)
x_max, x2, xd = int(sqrt((end-1)/4.0)), 0, 4
for xd in xrange(4, 8*x_max + 2, 8):
x2 += xd
y_max = int(sqrt(end-x2))
n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
if not (n & 1):
n -= n_diff
n_diff -= 2
for d in xrange((n_diff - 1) << 1, -1, -8):
m = n % 12
if m == 1 or m == 5:
m = n >> 1
sieve[m] = not sieve[m]
n -= d
x_max, x2, xd = int(sqrt((end-1) / 3.0)), 0, 3
for xd in xrange(3, 6 * x_max + 2, 6):
x2 += xd
y_max = int(sqrt(end-x2))
n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
if not(n & 1):
n -= n_diff
n_diff -= 2
for d in xrange((n_diff - 1) << 1, -1, -8):
if n % 12 == 7:
m = n >> 1
sieve[m] = not sieve[m]
n -= d
x_max, y_min, x2, xd = int((2 + sqrt(4-8*(1-end)))/4), -1, 0, 3
for x in xrange(1, x_max + 1):
x2 += xd
xd += 6
if x2 >= end: y_min = (((int(ceil(sqrt(x2 - end))) - 1) << 1) - 2) << 1
n, n_diff = ((x*x + x) << 1) - 1, (((x-1) << 1) - 2) << 1
for d in xrange(n_diff, y_min, -8):
if n % 12 == 11:
m = n >> 1
sieve[m] = not sieve[m]
n += d
primes = [2, 3]
if end <= 3:
return primes[:max(0,end-2)]
for n in xrange(5 >> 1, (int(sqrt(end))+1) >> 1):
if sieve[n]:
primes.append((n << 1) + 1)
aux = (n << 1) + 1
aux *= aux
for k in xrange(aux, end, 2 * aux):
sieve[k >> 1] = False
s = int(sqrt(end)) + 1
if s % 2 == 0:
s += 1
primes.extend([i for i in xrange(s, end, 2) if sieve[i >> 1]])
return primes
def ambi_sieve_plain(n):
s = range(3, n, 2)
for m in xrange(3, int(n**0.5)+1, 2):
if s[(m-3)/2]:
for t in xrange((m*m-3)/2,(n>>1)-1,m):
s[t]=0
return [2]+[t for t in s if t>0]
def sundaram3(max_n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/2073279#2073279
numbers = range(3, max_n+1, 2)
half = (max_n)//2
initial = 4
for step in xrange(3, max_n+1, 2):
for i in xrange(initial, half, step):
numbers[i-1] = 0
initial += 2*(step+1)
if initial > half:
return [2] + filter(None, numbers)
################################################################################
# Using Numpy:
def ambi_sieve(n):
# http://tommih.blogspot.com/2009/04/fast-prime-number-generator.html
s = np.arange(3, n, 2)
for m in xrange(3, int(n ** 0.5)+1, 2):
if s[(m-3)/2]:
s[(m*m-3)/2::m]=0
return np.r_[2, s[s>0]]
def primesfrom3to(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Returns a array of primes, p < n """
assert n>=2
sieve = np.ones(n/2, dtype=np.bool)
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i/2]:
sieve[i*i/2::i] = False
return np.r_[2, 2*np.nonzero(sieve)[0][1::]+1]
def primesfrom2to(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Input n>=6, Returns a array of primes, 2 <= p < n """
sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
sieve[0] = False
for i in xrange(int(n**0.5)/3+1):
if sieve[i]:
k=3*i+1|1
sieve[ ((k*k)/3) ::2*k] = False
sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]
if __name__=='__main__':
import itertools
import sys
def test(f1,f2,num):
print('Testing {f1} and {f2} return same results'.format(
f1=f1.func_name,
f2=f2.func_name))
if not all([a==b for a,b in itertools.izip_longest(f1(num),f2(num))]):
sys.exit("Error: %s(%s) != %s(%s)"%(f1.func_name,num,f2.func_name,num))
n=1000000
test(sieveOfAtkin,sieveOfEratosthenes,n)
test(sieveOfAtkin,ambi_sieve,n)
test(sieveOfAtkin,ambi_sieve_plain,n)
test(sieveOfAtkin,sundaram3,n)
test(sieveOfAtkin,sieve_wheel_30,n)
test(sieveOfAtkin,primesfrom3to,n)
test(sieveOfAtkin,primesfrom2to,n)
test(sieveOfAtkin,rwh_primes,n)
test(sieveOfAtkin,rwh_primes1,n)
test(sieveOfAtkin,rwh_primes2,n)
运行脚本测试所有实现均给出相同结果。
评论
如果您对非纯Python代码感兴趣，则应查看gmpy －通过mpz类型的next_prime方法，它对质数有很好的支持。
– Alex Martelli
2010-1-15的1:41
如果您使用的是pypy，那么这些基准（psyco基准）似乎还算不错。出乎意料的是，我发现sieveOfEratosthenes和ambi_sieve_plain是使用pypy最快的。这是我发现的非numpy的gist.github.com/5bf466bb1ee9e5726a52
– Ehsan Kia
2014年4月28日在1:53
如果有人想知道这里的功能与纯python而没有psyco或pypy的Wikibook的PG7.8相比如何：n = 1000000：PG7.8：每个循环4.93 s； rwh_primes1：每个循环69毫秒； rwh_primes2：每个循环57.1毫秒
–很棒
2015年7月9日在13:19
既然psyco已死并且PyPy已取代它，您可以用PyPy更新它吗？
–noɥʇʎԀʎzɐɹƆ
16年12月23日在22:20
如果可以针对python3更新这些功能和计时，那将是很棒的。
–cs95
18年11月12日下午5:22
#2 楼
更快，更明智的纯Python代码：def primes(n):
""" Returns a list of primes < n """
sieve = [True] * n
for i in range(3,int(n**0.5)+1,2):
if sieve[i]:
sieve[i*i::2*i]=[False]*((n-i*i-1)//(2*i)+1)
return [2] + [i for i in range(3,n,2) if sieve[i]]
或以半筛开始
def primes1(n):
""" Returns a list of primes < n """
sieve = [True] * (n//2)
for i in range(3,int(n**0.5)+1,2):
if sieve[i//2]:
sieve[i*i//2::i] = [False] * ((n-i*i-1)//(2*i)+1)
return [2] + [2*i+1 for i in range(1,n//2) if sieve[i]]
更快以及更多基于内存的numpy代码：
import numpy
def primesfrom3to(n):
""" Returns a array of primes, 3 <= p < n """
sieve = numpy.ones(n//2, dtype=numpy.bool)
for i in range(3,int(n**0.5)+1,2):
if sieve[i//2]:
sieve[i*i//2::i] = False
return 2*numpy.nonzero(sieve)[0][1::]+1
从三分之一的筛子开始的更快的变化：
import numpy
def primesfrom2to(n):
""" Input n>=6, Returns a array of primes, 2 <= p < n """
sieve = numpy.ones(n//3 + (n%6==2), dtype=numpy.bool)
for i in range(1,int(n**0.5)//3+1):
if sieve[i]:
k=3*i+1|1
sieve[ k*k//3 ::2*k] = False
sieve[k*(k-2*(i&1)+4)//3::2*k] = False
return numpy.r_[2,3,((3*numpy.nonzero(sieve)[0][1:]+1)|1)]
上述代码的（很难编码的）纯python版本是：
def primes2(n):
""" Input n>=6, Returns a list of primes, 2 <= p < n """
n, correction = n-n%6+6, 2-(n%6>1)
sieve = [True] * (n//3)
for i in range(1,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-2*(i&1)+4)//3::2*k] = [False] * ((n//6-k*(k-2*(i&1)+4)//6-1)//k+1)
return [2,3] + [3*i+1|1 for i in range(1,n//3-correction) if sieve[i]]
不幸的是，纯python不采用更简单和更快的numpy分配方式，并像在
len()中那样在循环内调用[False]*len(sieve[((k*k)//3)::2*k])太慢了。所以我不得不即兴改正输入（避免更多的数学运算），并做一些极端的（令人痛苦的）数学魔术。 我个人认为numpy（已被广泛使用）不是Python标准库的一部分是可耻的，并且Python开发人员似乎完全忽略了语法和速度方面的改进。 br />
评论
Numpy现在与Python 3兼容。事实上，它不在标准库中，因此它们可以拥有自己的发布周期。
–亚当
13年3月24日在23:17
只是将二进制值存储在数组中，我建议使用位数组-如此处所用（用于最简单的主筛；此处不是竞争者！）stackoverflow.com/questions/31120986/…
– hiro主角
15年7月17日在19:50
在使用primesfrom2to（）方法进行转换时，分隔是否应该在方括号内？
–355durch113
15-10-16在22:52
对于与python 3兼容的纯python版本，请点击以下链接：stackoverflow.com/a/33356284/2482582
– Moebius
16年4月8日在20:15
神圣的屁股这个吸盘很快。
–他
19年8月2日在7:52
#3 楼
这里有一个来自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()))
在pri.py中使用此代码在shell提示符下进行测量（我喜欢这样做），我观察到：
$ 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
所以看起来Cookbook解决方案的速度快了两倍。
评论
@jbochi，不客气-但请务必注意该URL，包括积分：我们花了十个人才将代码精炼到这一点，其中包括Tim Peters和Raymond Hettinger这样的Python性能专家（我写了自从我编辑印刷的Cookbook以来，食谱的最终文本，但是在编码方面，我的贡献与其他人差不多。）-最后，它确实是微妙且经过微调的代码，这不足为奇！-）
– Alex Martelli
2010年1月14日23:59
@Alex：知道您的代码“仅”是我的两倍，这让我感到非常自豪。 :) URL也非常有趣。再次感谢。
– jbochi
2010年1月15日，0：13
只需稍作更改，它就可以变得更快：请参见stackoverflow.com/questions/2211990/…
–tzot
2010-09-26在3:02
...可以通过延迟添加额外的〜1.2x-1.3x的速度，从O（n）到O（sqrt（n））显着减少内存占用量以及提高经验时间复杂度来使其更快。将prim填入dict，直到在输入中看到其平方为止。在这里测试。
–尼斯
2012年8月2日在22:28
#4 楼
我想使用Sundaram的Sieve打破了纯Python的记录：def sundaram3(max_n):
numbers = range(3, max_n+1, 2)
half = (max_n)//2
initial = 4
for step in xrange(3, max_n+1, 2):
for i in xrange(initial, half, step):
numbers[i-1] = 0
initial += 2*(step+1)
if initial > half:
return [2] + filter(None, numbers)
比较：
C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.get_primes_erat(1000000)"
10 loops, best of 3: 710 msec per loop
C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.daniel_sieve_2(1000000)"
10 loops, best of 3: 435 msec per loop
C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.sundaram3(1000000)"
10 loops, best of 3: 327 msec per loop
评论
我设法通过在函数顶部添加“ zero = 0”，然后将过滤器中的lambda替换为“ zero .__ sub__”，将函数的速度提高了约20％。不是世界上最漂亮的代码，但是速度更快：)
– truppo
2010年1月20日10:34
@truppo：感谢您的评论！我只是意识到传递None而不是原始函数是可行的，它甚至比零还快。
– jbochi
2010年1月20日，11：08
您是否知道如果您通过sundaram3（9），它将返回[2，3，5，7，9]？似乎用很多-也许所有-奇数（即使它们不是素数）也可以做到这一点
– wrhall
13年9月21日在23:57
它有一个问题：sundaram3（7071）包含7071，但不是素数
–bigOther
16年1月7日在7:08
#5 楼
该算法速度很快，但是存在严重缺陷：>>> sorted(get_primes(530))
[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, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163,
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 527, 529]
>>> 17*31
527
>>> 23*23
529
您假设
numbers.pop()将返回集合中的最小数字，但是完全不能保证。集是无序的，并且pop()删除并返回任意元素，因此不能用于从剩余数字中选择下一个素数。#6 楼
对于具有足够大N的真正最快的解决方案，是下载预先计算的素数列表，将其存储为元组，然后执行以下操作：for pos,i in enumerate(primes):
if i > N:
print primes[:pos]
如果仅
N > primes[-1]然后计算更多素数并将新列表保存在您的代码中，因此下次同样快。始终在框外思考。
评论
为了公平起见，您必须计算下载，解压缩和格式化素数的时间，然后将其与使用算法生成素数的时间进行比较-这些算法中的任何一种都可以轻松地将结果写入文件以供以后使用采用。我认为在这种情况下，如果有足够的内存来实际计算所有小于982,451,653的素数，那么numpy解决方案仍然会更快。
–丹尼尔G
2010-1-15在8:49
@Daniel正确。但是，您所拥有的商店仍然存在，并在需要时继续经营...
–金威
2010年1月15日上午10:11
@Daniel G我认为下载时间无关紧要。生成数字不是真的吗，所以您需要考虑用于创建要下载的列表的算法。给定时间O（n），任何时间复杂度都会忽略一次文件传输。
–罗斯
2014年10月22日22:56
UTM质数页面的常见问题解答建议，计算小质数比从磁盘读取小质数更快（问题是小意味什么）。
–蝙蝠侠
17年5月23日在0:59
#7 楼
如果您不想重新发明轮子，可以安装符号数学库sympy（是的，它与Python 3兼容）pip install sympy
并使用primerange函数
from sympy import sieve
primes = list(sieve.primerange(1, 10**6))
评论
我注意到这会打印整个列表，而社区Wiki答案primesfrom2to（10000）返回[2 3 5 ... 9949 9967 9973]。这会缩短NumPy nd.array吗？
– lukejanicke
8月4日19:15
#8 楼
如果您接受itertools但不接受numpy，则这是针对Python 3的rwh_primes2的改编版，在我的计算机上运行速度约为后者的两倍。唯一的实质性更改是使用字节数组而不是布尔值列表，并使用compress而不是列表推导来构建最终列表。 （如果可以的话，我将其添加为类似moarningsun的注释。）import itertools
izip = itertools.zip_longest
chain = itertools.chain.from_iterable
compress = itertools.compress
def rwh_primes2_python3(n):
""" Input n>=6, Returns a list of primes, 2 <= p < n """
zero = bytearray([False])
size = n//3 + (n % 6 == 2)
sieve = bytearray([True]) * size
sieve[0] = False
for i in range(int(n**0.5)//3+1):
if sieve[i]:
k=3*i+1|1
start = (k*k+4*k-2*k*(i&1))//3
sieve[(k*k)//3::2*k]=zero*((size - (k*k)//3 - 1) // (2 * k) + 1)
sieve[ start ::2*k]=zero*((size - start - 1) // (2 * k) + 1)
ans = [2,3]
poss = chain(izip(*[range(i, n, 6) for i in (1,5)]))
ans.extend(compress(poss, sieve))
return ans
比较：
>>> timeit.timeit('primes.rwh_primes2(10**6)', setup='import primes', number=1)
0.0652179726976101
>>> timeit.timeit('primes.rwh_primes2_python3(10**6)', setup='import primes', number=1)
0.03267321276325674
和
>>> timeit.timeit('primes.rwh_primes2(10**8)', setup='import primes', number=1)
6.394284538007014
>>> timeit.timeit('primes.rwh_primes2_python3(10**8)', setup='import primes', number=1)
3.833829450302801
#9 楼
编写自己的主要发现代码很有启发性，但是手头有一个快速可靠的库也很有用。我围绕C ++库primesieve写了一个包装，将其命名为primesieve-python试试看
pip install primesieve import primesieve
primes = primesieve.generate_primes(10**8)
我很想知道比较的速度。
评论
这不完全是OP的命令，但我看不出为什么要投票。与其他一些外部模块不同，这是一个2.8秒的解决方案。我在源中注意到它是线程化的，是否对其扩展性进行了任何测试？
– ljetibo
15年7月14日在7:24
@ljetibo欢呼。瓶颈似乎是将C ++向量复制到Python列表，因此count_primes函数比generate_primes快得多
–上校恐慌
15年7月14日在8:19
在我的计算机上，它可以舒适地生成最多为1e8的素数（它为1e9提供MemoryError），并可以计算最多为1e10的素数。上面的@HappyLeapSecond比较了1e6的算法
–上校恐慌
15年7月14日在8:30
#10 楼
这是最快的功能之一的两个更新版本（纯Python 3.6）。from itertools import compress
def rwh_primes1v1(n):
""" Returns a list of primes < n for n > 2 """
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)
return [2,*compress(range(3,n,2), sieve[1:])]
def rwh_primes1v2(n):
""" Returns a list of primes < n for n > 2 """
sieve = bytearray([True]) * (n//2+1)
for i in range(1,int(n**0.5)//2+1):
if sieve[i]:
sieve[2*i*(i+1)::2*i+1] = bytearray((n//2-2*i*(i+1))//(2*i+1)+1)
return [2,*compress(range(3,n,2), sieve[1:])]
评论
在Python 3中，我使用了这个函数stackoverflow.com/a/3035188/7799269，但用/替换了/，并用range替换了xrange，它们看起来比这快得多。
–samerivertwice
18年9月11日在21:20
#11 楼
根据N <9,080,191import sys
import random
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 xrange(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):
for a in [2, 3, 37, 73]:
if not miller_rabin_pass(a, n):
return False
return True
n = int(sys.argv[1])
primes = [2]
for p in range(3,n,2):
if miller_rabin(p):
primes.append(p)
print len(primes)
的确定性实现Miller-Rabin素数检验的确定性实现（根据Wikipedia上的文章（http://en.wikipedia.org/ wiki / Miller–Rabin_primality_test）测试N <9,080,191的a = 2,3,37和73足以决定N是否为复合。
我从概率实现中改编了源代码可在此处找到原始Miller-Rabin的测试：http://en.literateprograms.org/Miller-Rabin_primality_test_(Python）
评论
感谢Miller-Rabin素数测试，但是该代码实际上速度较慢，并且无法提供正确的结果。 37是质数，未通过测试。
– jbochi
2010年1月21日在10:07
我想37是特殊情况之一，我不好。我对确定性版本抱有希望:)
– Ruggiero Spearman
2010年1月21日，11：07
Rabin Miller没有任何特殊情况。
–错误
13年11月26日在23:37
您看错了这篇文章。它是31，而不是37。这就是为什么您的实现失败的原因。
–登录
13年12月24日在22:36
#12 楼
如果您可以控制N，则列出所有素数的最快方法是预先计算它们。说真的预计算是一种被忽略的优化方法。评论
或从此处primes.utm.edu/lists/small/millions下载它们，但其想法是测试python的限制，并查看优化是否产生了漂亮的代码。
– jbochi
2010年1月21日，10：14
#13 楼
这是我通常用于在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。 br />
感谢您发布此问题。今天我真的学到了很多东西。
#14 楼
对于最快的代码，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
#15 楼
使用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半筛媲美。
评论
最高= 10 ** 6：primesfrom2to（）-7毫秒; prime6（）-12毫秒ideone.com/oDg2Y
– jfs
10-9-4'1:57
#16 楼
全部都经过编写和测试。因此，无需重新设计轮子。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毫秒的最佳解决方案。
比较一下！！
#17 楼
我测试了一些unutbu的函数，并用数以百万计的数字进行了计算赢家是使用numpy库的函数，
注意：进行内存利用率测试也很有趣:)
示例代码
我的github存储库中的完整代码
#!/usr/bin/env python
import lib
import timeit
import sys
import math
import datetime
import prettyplotlib as ppl
import numpy as np
import matplotlib.pyplot as plt
from prettyplotlib import brewer2mpl
primenumbers_gen = [
'sieveOfEratosthenes',
'ambi_sieve',
'ambi_sieve_plain',
'sundaram3',
'sieve_wheel_30',
'primesfrom3to',
'primesfrom2to',
'rwh_primes',
'rwh_primes1',
'rwh_primes2',
]
def human_format(num):
# https://stackoverflow.com/questions/579310/formatting-long-numbers-as-strings-in-python?answertab=active#tab-top
magnitude = 0
while abs(num) >= 1000:
magnitude += 1
num /= 1000.0
# add more suffixes if you need them
return '%.2f%s' % (num, ['', 'K', 'M', 'G', 'T', 'P'][magnitude])
if __name__=='__main__':
# Vars
n = 10000000 # number itereration generator
nbcol = 5 # For decompose prime number generator
nb_benchloop = 3 # Eliminate false positive value during the test (bench average time)
datetimeformat = '%Y-%m-%d %H:%M:%S.%f'
config = 'from __main__ import n; import lib'
primenumbers_gen = {
'sieveOfEratosthenes': {'color': 'b'},
'ambi_sieve': {'color': 'b'},
'ambi_sieve_plain': {'color': 'b'},
'sundaram3': {'color': 'b'},
'sieve_wheel_30': {'color': 'b'},
# # # 'primesfrom2to': {'color': 'b'},
'primesfrom3to': {'color': 'b'},
# 'rwh_primes': {'color': 'b'},
# 'rwh_primes1': {'color': 'b'},
'rwh_primes2': {'color': 'b'},
}
# Get n in command line
if len(sys.argv)>1:
n = int(sys.argv[1])
step = int(math.ceil(n / float(nbcol)))
nbs = np.array([i * step for i in range(1, int(nbcol) + 1)])
set2 = brewer2mpl.get_map('Paired', 'qualitative', 12).mpl_colors
print datetime.datetime.now().strftime(datetimeformat)
print("Compute prime number to %(n)s" % locals())
print("")
results = dict()
for pgen in primenumbers_gen:
results[pgen] = dict()
benchtimes = list()
for n in nbs:
t = timeit.Timer("lib.%(pgen)s(n)" % locals(), setup=config)
execute_times = t.repeat(repeat=nb_benchloop,number=1)
benchtime = np.mean(execute_times)
benchtimes.append(benchtime)
results[pgen] = {'benchtimes':np.array(benchtimes)}
fig, ax = plt.subplots(1)
plt.ylabel('Computation time (in second)')
plt.xlabel('Numbers computed')
i = 0
for pgen in primenumbers_gen:
bench = results[pgen]['benchtimes']
avgs = np.divide(bench,nbs)
avg = np.average(bench, weights=nbs)
# Compute linear regression
A = np.vstack([nbs, np.ones(len(nbs))]).T
a, b = np.linalg.lstsq(A, nbs*avgs)[0]
# Plot
i += 1
#label="%(pgen)s" % locals()
#ppl.plot(nbs, nbs*avgs, label=label, lw=1, linestyle='--', color=set2[i % 12])
label="%(pgen)s avg" % locals()
ppl.plot(nbs, a * nbs + b, label=label, lw=2, color=set2[i % 12])
print datetime.datetime.now().strftime(datetimeformat)
ppl.legend(ax, loc='upper left', ncol=4)
# Change x axis label
ax.get_xaxis().get_major_formatter().set_scientific(False)
fig.canvas.draw()
labels = [human_format(int(item.get_text())) for item in ax.get_xticklabels()]
ax.set_xticklabels(labels)
ax = plt.gca()
plt.show()
评论
为了比较算法性能，最好以对数-对数比例绘制。
–尼斯
19年8月26日在9:35
#18 楼
对于Python 3def 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]]
#19 楼
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
我优化了Eratosthenes筛，以提高速度和内存。
基准
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
#20 楼
第一次使用python，因此我在其中使用的某些方法似乎有点麻烦。我只是将我的c ++代码直接转换为python，这就是我所拥有的（尽管python中有一点点slowww）#!/usr/bin/env python
import time
def GetPrimes(n):
Sieve = [1 for x in xrange(n)]
Done = False
w = 3
while not Done:
for q in xrange (3, n, 2):
Prod = w*q
if Prod < n:
Sieve[Prod] = 0
else:
break
if w > (n/2):
Done = True
w += 2
return Sieve
start = time.clock()
d = 10000000
Primes = GetPrimes(d)
count = 1 #This is for 2
for x in xrange (3, d, 2):
if Primes[x]:
count+=1
elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"
pythonw Primes.py
在12.799119秒内发现664579个质数！
#!/usr/bin/env python
import time
def GetPrimes2(n):
Sieve = [1 for x in xrange(n)]
for q in xrange (3, n, 2):
k = q
for y in xrange(k*3, n, k*2):
Sieve[y] = 0
return Sieve
start = time.clock()
d = 10000000
Primes = GetPrimes2(d)
count = 1 #This is for 2
for x in xrange (3, d, 2):
if Primes[x]:
count+=1
elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"
pythonw Primes2.py
在10.230172秒内找到664579素数！
#!/usr/bin/env python
import time
def GetPrimes3(n):
Sieve = [1 for x in xrange(n)]
for q in xrange (3, n, 2):
k = q
for y in xrange(k*k, n, k << 1):
Sieve[y] = 0
return Sieve
start = time.clock()
d = 10000000
Primes = GetPrimes3(d)
count = 1 #This is for 2
for x in xrange (3, d, 2):
if Primes[x]:
count+=1
elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"
python Primes2.py
在7.113776中找到664579素数。秒！
#21 楼
我知道比赛已经结束了几年。 …尽管如此，我还是建议您使用纯python素数筛，该方法基于在处理正数筛时通过使用适当的步骤来省略2、3和5的倍数。但是，对于N <10 ^ 9而言，它实际上要比@Robert William Hanks高级解决方案rwh_primes2和rwh_primes1慢。通过使用高于1.5 * 10 ^ 8的ctypes.c_ushort筛子数组，可以以某种方式适应内存限制。
10 ^ 6
$ python -mtimeit -s“ import primeSieveSpeedComp “” primeSieveSpeedComp.primeSieveSeq（1000000）“
10个循环，每个循环最好3：46.7毫秒
要比较：$ python -mtimeit -s” import primeSieveSpeedComp“
“ primeSieveSpeedComp.rwh_primes1（1000000）” 10个循环，最好为3：43.2
每个循环毫秒[br />比较：$ python -m timeit -s“ import primeSieveSpeedComp”
“ primeSieveSpeedComp .rwh_primes2（1000000）“ 10个循环，最好为3：34.5
每个循环毫秒（
10 ^ 7
$ python -mtimeit -s “ import primeSieveSpeedComp”“ primeSieveSpeedComp.primeSieveSeq（10000000）”
10个循环，每循环最好3：530毫秒
进行比较：$ python -mtimeit -s“ import primeSieveSpeedComp“
” primeSieveSpeedComp.rwh_primes1（10000000）“ 10个循环，最好为3：494
每个循环毫秒（
）进行比较：$ python -m timeit -s” import primeSieveSpeedComp“”
“ primeSieveSpeedComp.rwh_primes2（10000000）”最佳10循环，最好3：375
每个循环msec
10 ^ 8
$ python -mtimeit -s“ import primeSieveSpeedComp”“ primeSieveSpeedComp.primeSieveSeq（100000000）”
10个循环，最好是3：每个循环5.55秒
比较： $ python -mtimeit -s“ import primeSieveSpeedComp”
“ primeSieveSpeedComp.rwh_primes1（100000000）” 10个循环，最好是3：5.33
每个循环秒
比较：$ python -m timeit- s“ import primeSieveSpeedComp”
“ primeSieveSpeedComp.rwh_primes2（100000000）” 10个循环，最好为3：3.95
每个循环秒数
10 ^ 9
$ python -mtimeit -s“ import primeSieveSpeedComp”“ primeSieveSpeedComp.primeSieveSeq（1000000000）”
10个循环，最好是3：每个循环61.2秒
比较：$ python -mtimeit -n 3 -s“ import primeSieveSpeedComp”
“ primeSieveSpeedComp.rwh_primes1（1000000000）” 3循环，最好为3：每个循环97.8
秒
比较：$ python -m timeit -s“ import primeSieveSpeedComp”
“ primeSieveSpeedComp.rwh_primes2（1000000000）” 10个循环，效果最好：3：
每个循环41.9秒
您可以将以下代码复制到ubuntus primeSieveSpeedComp中，以查看此测试。
def primeSieveSeq(MAX_Int):
if MAX_Int > 5*10**8:
import ctypes
int16Array = ctypes.c_ushort * (MAX_Int >> 1)
sieve = int16Array()
#print 'uses ctypes "unsigned short int Array"'
else:
sieve = (MAX_Int >> 1) * [False]
#print 'uses python list() of long long int'
if MAX_Int < 10**8:
sieve[4::3] = [True]*((MAX_Int - 8)/6+1)
sieve[12::5] = [True]*((MAX_Int - 24)/10+1)
r = [2, 3, 5]
n = 0
for i in xrange(int(MAX_Int**0.5)/30+1):
n += 3
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 2
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 1
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 2
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 1
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 2
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 3
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 1
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
if MAX_Int < 10**8:
return [2, 3, 5]+[(p << 1) + 1 for p in [n for n in xrange(3, MAX_Int >> 1) if not sieve[n]]]
n = n >> 1
try:
for i in xrange((MAX_Int-2*n)/30 + 1):
n += 3
if not sieve[n]:
r.append((n << 1) + 1)
n += 2
if not sieve[n]:
r.append((n << 1) + 1)
n += 1
if not sieve[n]:
r.append((n << 1) + 1)
n += 2
if not sieve[n]:
r.append((n << 1) + 1)
n += 1
if not sieve[n]:
r.append((n << 1) + 1)
n += 2
if not sieve[n]:
r.append((n << 1) + 1)
n += 3
if not sieve[n]:
r.append((n << 1) + 1)
n += 1
if not sieve[n]:
r.append((n << 1) + 1)
except:
pass
return r
评论
要可视化您的测试结果，以对数-对数比例绘制它们，以查看和比较增长的经验顺序。
–尼斯
16-3-3在12:18
@感谢您的输入，下次我需要这种比较时，请记住这一点
– ABri
16年3月12日在20:04
#22 楼
我已经更新了许多Python 3代码，并将其扔到了perfplot（我的一个项目）上，看看哪个实际上最快。事实证明，对于较大的n，primesfrom{2,3}to占了上风：
代码重现该图：
#23 楼
这是Eratosthenes筛子的一个小版本，具有很好的复杂性（比对长度为n的数组排序要低）和矢量化。与@unutbu相比，此速度与包含46个微秒的包装的速度一样快，可以找到所有百万以下的素数。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
#24 楼
我的猜测是，所有方法中最快的就是对代码中的质数进行硬编码。所以为什么不编写一个慢速脚本来生成另一个包含所有数字的源文件，然后当您运行实际程序时，请导入该源文件。
当然，只有在编译时知道N的上限时，此方法才有效，但是（几乎）所有项目Euler问题都是这种情况。
PS：我可能错了，尽管使用硬连线素数解析源代码比首先计算它们要慢，但是据我所知Python运行于已编译了
.pyc文件，因此在这种情况下，读取所有素数最多为N的二进制数组应该非常快。#25 楼
很抱歉打扰，但是erat2（）在算法中存在严重缺陷。在搜索下一个组合时，我们只需要测试奇数。
q，p都是奇数；那么q + p是偶数，不需要测试，但是q + 2 * p总是奇数。这消除了while循环条件下的“ if even”测试，并节省了大约30％的运行时间。
我们正在这样做：而不是优雅的'D.pop（q，None） 'get and delete方法使用'if D in q：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
评论
有关将素数推迟添加到dict的信息（直到在输入中看到素数的平方），请参见stackoverflow.com/a/10733621/849891。
–尼斯
2012年11月10日21:12
#26 楼
到目前为止，我尝试过的最快方法是基于Python食谱erat2函数：import itertools as it
def erat2a( ):
D = { }
yield 2
for q in it.islice(it.count(3), 0, None, 2):
p = D.pop(q, None)
if p is None:
D[q*q] = q
yield q
else:
x = q + 2*p
while x in D:
x += 2*p
D[x] = p
请参见此答案以获取加速说明。
#27 楼
我可能聚会晚了，但是必须为此添加我自己的代码。它占用的空间大约为n / 2，因为我们不需要存储偶数，而且我还使用了bitarray python模块，从而进一步大大减少了内存消耗，并使所有素数的计算量高达1,000,000,000from 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
这是在64位2.4GHZ MAC OSX 10.8.3上运行的
评论
为未知机器发布一个时序不会说明任何事情。此处接受的答案是“没有psyco，对于n = 1000000，rwh_primes2是最快的”。因此，如果您要在同一台机器上同时提供该代码和您的代码的时间，并且同时提供2，4，10百万，那么它将提供更多信息。
–尼斯
13年4月17日在7:32
-1，此代码取决于用C实现的位数组的特殊功能，这就是为什么代码如此之快的原因，因为大部分工作是在分片分配的本机代码中完成的。位数组程序包打破了可变序列正确切片的标准定义（在一定范围内建立索引），因为它允许为切片的所有元素分配一个布尔值0/1或True / False，而纯Python的标准行为似乎是不允许这样做，只允许赋值为0，在这种情况下，它将其视为序列/数组中所有切片元素的del。
–GordonBGood
13年8月19日在16:46
续：如果要比较调用非标准的本机代码，我们不妨编写一个基于C代码的“ fastprimes”序列生成器程序包，例如Kim Walisch的primesieve，并生成40亿加32中的所有素数只需调用序列生成器，即可在几秒钟之内获得-位数字范围。由于链接的代码基于分段的Eratosthenes筛网，因此这也几乎不使用内存，因此仅使用了几十KB的RAM，并且如果生成了序列，则不需要列表存储。
–GordonBGood
13年8月19日在17:54
#28 楼
随着时间的推移，我收集了几个质数筛。我计算机上最快的计算机是：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
#29 楼
我对这个问题的回答很慢，但这似乎很有趣。我正在使用numpy，这可能会作弊，我怀疑这种方法是最快的，但应该清楚。它筛选仅引用其索引的布尔数组，并从所有True值的索引中引出质数。无需取模。import numpy as np
def ajs_primes3a(upto):
mat = np.ones((upto), dtype=bool)
mat[0] = False
mat[1] = False
mat[4::2] = False
for idx in range(3, int(upto ** 0.5)+1, 2):
mat[idx*2::idx] = False
return np.where(mat == True)[0]
评论
这是不正确的，例如ajs_primes3a（10）-> array（[2，3，5，7，9]）。 9不是素数
– jfs
15年2月14日在21:06
您发现了一个我还没有的尖端案例-做得好！问题出在'对于range（3，int（最大** 0.5），2）中的idx：'应该是'对于range（3，int（最大** 0.5）+ 1，2）中的idx：'。谢谢，但是现在可以使用。
–艾伦·詹姆斯·萨尔莫尼（Alan James Salmoni）
2015年2月15日在22:14
原因是idx循环上升到“最高** 05”，对于15或以下的情况，从16开始，它可以正常工作。这是我未测试过的一组边缘案例。加1表示它应适用于所有数字。
–艾伦·詹姆斯·萨尔莫尼（Alan James Salmoni）
2015年2月15日在22:18
现在似乎可以正常工作了。在返回数组的基于numpy的解决方案中，它是最慢的。注意：Eratosthenes实现的真正筛子不使用模数-无需提及。您可以使用mat [idx * idx :: idx]代替mat [idx * 2 :: idx]。和np.nonzero（mat）[0]而不是np.where（mat == True）[0]。
– jfs
2015年2月15日在23:53
谢谢JF。我针对prime6（）进行了测试，当prime6（）接手时，得到的结果更快地达到了（IIRC）约250k。 primesfrom2to（）更快。在最长20m的时间内，ajs_primes3a（）占用0.034744977951ms，prime6（）占用0.0222899913788ms，primesfrom2to（）占用0.0104751586914ms（同一台机器，相同的负载，最好是10个计时）。说实话，这比我想象的要好！
–艾伦·詹姆斯·萨尔莫尼（Alan James Salmoni）
15年2月16日在15:54
#30 楼
这是一种使用python的列表推导来生成质数（但不是最有效的）的有趣技术：noprimes = [j for i in range(2, 8) for j in range(i*2, 50, i)]
primes = [x for x in range(2, 50) if x not in noprimes]
您可以在此处找到示例和一些说明
评论
如果像数字一样声明的数字= set（range（n，2，-2）），则所讨论的代码段将更快。但是不能击败圣达拉姆。感谢您的提问。如果答案中可以有Python 3版本的函数，那将是很好的。
当然，有一个库可以执行此操作，因此我们不必滚动我们自己的> xkcd，它保证Python就像导入反重力一样简单。是否没有需要“素数”的东西？ Prime.take（10）（Ruby）？
@ColonelPanic就这样，我为Py3更新了github.com/jaredks/pyprimesieve，并添加到了PyPi中。它肯定比这些要快，但幅度不大-比最佳的numpy版本快约5倍。
@ColonelPanic：我认为编辑旧答案以指出它们已经老化是适当的，因为这使它成为更有用的资源。如果“已接受”的答案不再是最佳答案，则可以在问题中添加注释，并在2015年进行更新以指出人们当前的最佳方法。