列出 N 以下所有素数的最快方法 [英] Fastest way to list all primes below N
本文介绍了列出 N 以下所有素数的最快方法的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
这是我能想到的最好的算法.
def get_primes(n):数字 = 设置(范围(n，1，-1))素数 = []而数字:p = numbers.pop()素数.append(p)numbers.difference_update(set(range(p*2,n+1,p)))返回素数>>>timeit.Timer(stmt='get_primes.get_primes(1000000)', setup='import get_primes').timeit(1)1.1499958793645562可以做得更快吗?
这段代码有一个缺陷:由于numbers 是一个无序集合，所以不能保证numbers.pop() 会从集合中删除最小的数字.尽管如此，它对某些输入数字有效(至少对我而言):
解决方案
警告: timeit 结果可能因硬件或版本的 Python.
下面是一个脚本，它比较了许多实现:
- ambi_sieve_plain，
- rwh_primes,
- rwh_primes1,
- rwh_primes2,
- sieveOfAtkin,
- sieveOfEratosthenes,
- sundaram3,
- sieve_wheel_30,
- ambi_sieve(需要 numpy)
- primesfrom3to(需要麻木)
- primesfrom2to(需要麻木)
非常感谢 stephan 让我注意到了筛轮_30.归功于 罗伯特·威廉·汉克斯 对于 primesfrom2to、primesfrom3to、rwh_primes、rwh_primes1 和 rwh_primes2.
在测试的纯 Python 方法中，使用 psyco，对于 n=1000000，rwh_primes1 是经过测试的最快的.
+---------------------+-------+|方法 |女士 |+---------------------+-------+|rwh_primes1 |43.0 ||筛子阿特金 |46.4 ||rwh_primes |57.4 ||筛轮_30 |63.0 ||rwh_primes2 |67.8 ||埃拉托色尼筛 |147.0 ||ambi_sieve_plain |152.0 ||sundaram3 |194.0 |+---------------------+-------+在测试的纯 Python 方法中，没有 psyco，对于 n=1000000，rwh_primes2 是最快的.
+---------------------+-------+|方法 |女士 |+---------------------+-------+|rwh_primes2 |68.1 ||rwh_primes1 |93.7 ||rwh_primes |94.6 ||筛轮_30 |97.4 ||埃拉托色尼筛 |178.0 ||ambi_sieve_plain |286.0 ||筛子阿特金 |314.0 ||sundaram3 |416.0 |+---------------------+-------+在所有测试的方法中，允许 numpy，对于 n=1000000，primesfrom2to 是测试最快的.
+---------------------+-------+|方法 |女士 |+---------------------+-------+|素数从2到|15.9 ||素数从3到|18.4 ||ambi_sieve |29.3 |+---------------------+-------+使用以下命令测量时间:
python -mtimeit -s"import primes" "primes.{method}(1000000)"用 {method} 替换为每个方法名称.
primes.py:
#!/usr/bin/env python进口精神科；psyco.full()从数学导入 sqrt, ceil将 numpy 导入为 npdef rwh_primes(n):# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188""" 返回一个质数列表 < n """筛 = [真] * n对于 xrange(3,int(n**0.5)+1,2) 中的 i:如果筛[i]:筛[i*i::2*i]=[假]*((n-i*i-1)/(2*i)+1)返回 [2] + [i for i in xrange(3,n,2) 如果筛 [i]]def rwh_primes1(n):# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188""" 返回一个质数列表 < n """筛 = [真] * (n/2)对于 xrange(3,int(n**0.5)+1,2) 中的 i:如果筛[i/2]:筛[i*i/2::i] = [假] * ((n-i*i-1)/(2*i)+1)返回 [2] + [2*i+1 for i in xrange(1,n/2) 如果筛 [i]]def rwh_primes2(n):# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188""" 输入 n>=6，返回素数列表，2 <= p 1)n = {0:n,1:n-1,2:n+4,3:n+3,4:n+2,5:n+1}[n%6]筛 = [真] * (n/3)筛 [0] = 假对于 xrange(int(n**0.5)/3+1) 中的 i:如果筛[i]:k=3*i+1|1筛 [ ((k*k)/3) ::2*k]=[假]*((n/6-(k*k)/6-1)/k+1)筛[(k*k+4*k-2*k*(i&1))/3::2*k]=[假]*((n/6-(k*k+4*k-2*k*(i&1))/6-1)/k+1)返回 [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) 如果筛 [i]]定义筛轮_30(N):# http://zerovolt.com/?p=88''' 使用轮子标准 2*3*5 = 30 返回质数列表 <= N版权所有 2009 zerovolt.com此代码可免费用于非商业目的，在这种情况下，您可以将此评论作为我工作的功劳.如果您需要此代码用于商业目的，请发送电子邮件至: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)轮 = (2, 3, 5)常量 = 30如果 N＜2:返回 []如果 N <= 常量:位置 = 0而 __smallp[pos] <= N:位置 += 1返回列表(__smallp[:pos])# 制作偏移量列表偏移量 = (7, 11, 13, 17, 19, 23, 29, 1)# 准备清单p = [2, 3, 5]暗淡 = 2 + N//常量tk1 = [真] * 暗淡tk7 = [真] * 暗淡tk11 = [真] * 暗淡tk13 = [真] * 暗淡tk17 = [真] * 暗淡tk19 = [真] * 暗淡tk23 = [真] * 暗tk29 = [真] * 暗淡tk1[0] = 假# 帮助字典 d# d[a, b] = c ==>如果我想找到 (30*pos)+a 的最小有用倍数# 在 tkc 上，那么我需要由 [(30*pos)+a][(30*pos)+b] 的乘积给出的索引# 一般来说.如果 b＜a，我需要 [(30*pos)+a][(30*(pos+1))+b]d = {}对于偏移量中的 x:对于偏移量中的 y:res = (x*y) % 常量如果 res 为偏移量:d[(x, res)] = y# 另一个帮助字典:给出 tkx 调用 tmptk[x]tmptk = {1:tk1, 7:tk7, 11:tk11, 13:tk13, 17:tk17, 19:tk19, 23:tk23, 29:tk29}位置，素数，最后添加，停止 = 0, 0, 0, int(ceil(sqrt(N)))# 内部函数定义def del_mult(tk, start, step):对于 xrange(start, len(tk), step) 中的 k:tk[k] = 假# 内部函数定义结束cpos = const * pos而素数<停止:# 30k + 7如果 tk7[pos]:素数 = cpos + 7p.append(prime)最后添加 = 7在偏移量中关闭:tmp = d[(7, 关闭)]start = (pos + prime) if off == 7 else (prime * (const * (pos + 1 if tmp < 7 else 0) + tmp) )//constdel_mult(tmptk[off], 开始, 素数)# 30k + 11如果 tk11[pos]:素数 = cpos + 11p.append(prime)最后添加 = 11在偏移量中关闭:tmp = d[(11, 关闭)]start = (pos + prime) if off == 11 else (prime * (const * (pos + 1 if tmp < 11 else 0) + tmp) )//constdel_mult(tmptk[off], 开始, 素数)# 30k + 13如果 tk13[pos]:素数 = cpos + 13p.append(prime)最后添加 = 13在偏移量中关闭:tmp = d[(13, 关闭)]start = (pos + prime) if off == 13 else (prime * (const * (pos + 1 if tmp < 13 else 0) + tmp) )//constdel_mult(tmptk[off], 开始, 素数)# 30k + 17如果 tk17[pos]:素数 = cpos + 17p.append(prime)最后添加 = 17在偏移量中关闭:tmp = d[(17, 关闭)]start = (pos + prime) if off == 17 else (prime * (const * (pos + 1 if tmp < 17 else 0) + tmp) )//constdel_mult(tmptk[off], 开始, 素数)# 30k + 19如果 tk19[pos]:素数 = cpos + 19p.append(prime)最后添加 = 19在偏移量中关闭:tmp = d[(19, 关闭)]start = (pos + prime) if off == 19 else (prime * (const * (pos + 1 if tmp < 19 else 0) + tmp) )//constdel_mult(tmptk[off], 开始, 素数)# 30k + 23如果 tk23[pos]:素数 = cpos + 23p.append(prime)最后添加 = 23在偏移量中关闭:tmp = d[(23, 关闭)]start = (pos + prime) if off == 23 else (prime * (const * (pos + 1 if tmp < 23 else 0) + tmp) )//constdel_mult(tmptk[off], 开始, 素数)# 30k + 29如果 tk29[pos]:素数 = cpos + 29p.append(prime)最后添加 = 29在偏移量中关闭:tmp = d[(29, 关闭)]start = (pos + prime) if off == 29 else (prime * (const * (pos + 1 if tmp < 29 else 0) + tmp) )//constdel_mult(tmptk[off], 开始, 素数)# 现在我们回到顶部 tk1，所以我们需要将 pos 增加 1位置 += 1cpos = const * pos# 30k + 1如果 tk1[pos]:素数 = cpos + 1p.append(prime)最后添加 = 1在偏移量中关闭:tmp = d[(1, 关闭)]start = (pos + prime) if off == 1 else (prime * (const * pos + tmp) )//constdel_mult(tmptk[off], 开始, 素数)# 添加剩余素数的时间# 如果 last added == 1，删除最后一个元素并从 tk1 开始添加它们# 这样我们就不需要在最后一段时间内使用if"如果上次添加 == 1:p.pop()# 现在对所有其他可能的素数都完成当 pos <长度(tk1):cpos = const * pos如果 tk1[pos]: p.append(cpos + 1)如果 tk7[pos]: p.append(cpos + 7)如果 tk11[pos]: p.append(cpos + 11)如果 tk13[pos]: p.append(cpos + 13)如果 tk17[pos]: p.append(cpos + 17)如果 tk19[pos]: p.append(cpos + 19)如果 tk23[pos]: p.append(cpos + 23)如果 tk29[pos]: p.append(cpos + 29)位置 += 1# 删除超过如果存在pos = len(p) - 1而 p[pos] >否:位置 -= 1如果 pos > 1]])返回素数def ambi_sieve_plain(n):s = 范围(3, n, 2)对于 xrange(3, int(n**0.5)+1, 2) 中的 m:如果 s[(m-3)/2]:对于 xrange((m*m-3)/2,(n>>1)-1,m) 中的 t:s[t]=0返回 [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数字 = 范围(3，max_n+1，2)一半 = (max_n)//2初始 = 4对于 xrange(3, max_n+1, 2) 中的步骤:对于 xrange 中的 i(初始、一半、步长):数字[i-1] = 0初始 += 2*(步+1)如果初始 >一半:返回 [2] + 过滤器(无，数字)############################################################################# 使用 Numpy:def ambi_sieve(n):# http://tommih.blogspot.com/2009/04/fast-prime-number-generator.htmls = np.arange(3, n, 2)对于 xrange(3, int(n ** 0.5)+1, 2) 中的 m:如果 s[(m-3)/2]:s[(m*m-3)/2::m]=0返回 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""" 返回素数数组，p =2筛 = np.ones(n/2, dtype=np.bool)对于 xrange(3,int(n**0.5)+1,2) 中的 i:如果筛[i/2]:筛 [i*i/2::i] = 假返回 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""" 输入 n>=6, 返回素数数组, 2 <= p 运行脚本测试所有实现都给出相同的结果.
This is the best algorithm I could come up.
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
Can it be made even faster?
This code has a flaw: Since numbers is an unordered set, there is no guarantee that numbers.pop() will remove the lowest number from the set. Nevertheless, it works (at least for me) for some input numbers:
>>> 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
解决方案
Warning: timeit results may vary due to differences in hardware or
version of Python.
Below is a script which compares a number of implementations:
- ambi_sieve_plain,
- rwh_primes,
- rwh_primes1,
- rwh_primes2,
- sieveOfAtkin,
- sieveOfEratosthenes,
- sundaram3,
- sieve_wheel_30,
- ambi_sieve (requires numpy)
- primesfrom3to (requires numpy)
- primesfrom2to (requires numpy)
Many thanks to stephan for bringing sieve_wheel_30 to my attention. Credit goes to Robert William Hanks for primesfrom2to, primesfrom3to, rwh_primes, rwh_primes1, and rwh_primes2.
Of the plain Python methods tested, with psyco, for n=1000000, rwh_primes1 was the fastest tested.
+---------------------+-------+
| 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 |
+---------------------+-------+
Of the plain Python methods tested, without psyco, for n=1000000, rwh_primes2 was the fastest.
+---------------------+-------+
| 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 |
+---------------------+-------+
Of all the methods tested, allowing numpy, for n=1000000, primesfrom2to was the fastest tested.
+---------------------+-------+
| Method | ms |
+---------------------+-------+
| primesfrom2to | 15.9 |
| primesfrom3to | 18.4 |
| ambi_sieve | 29.3 |
+---------------------+-------+
Timings were measured using the command:
python -mtimeit -s"import primes" "primes.{method}(1000000)"
with {method} replaced by each of the method names.
primes.py:
#!/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)
Running the script tests that all implementations give the same result.
这篇关于列出 N 以下所有素数的最快方法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
