In [1]:
N = 100
table = [True] * N # 無罪推定, 先假定所有數字都是質數
table[0] = table[1] = False # 0 和 1 不是質數
for i in range(2,N):
if table[i]:
# 所有 i 的倍數都不是質數
for j in range(i*i, N, i):
table[j] = False
In [2]:
# 把結果列出來
for i in range(N):
print i, table[i]
測試一下速度
In [3]:
%%timeit
N = 10**7
table = [True] * N # 無罪推定, 先假定所有數字都是質數
table[0] = table[1] = False # 0 和 1 不是質數
for i in range(2,N):
if table[i]:
# 所有 i 的倍數都不是質數
for j in range(i*i, N, i):
table[j] = False
改進一些
In [12]:
%%timeit
N = 10**7
table = [False, True] * (N//2)
table[1:3] = [False, True]
for i in xrange(3,int(N**0.5)+1,2):
if table[i]:
i2 = i*i
table[i2::i] = [False]*((N-1-i2)//i+1)
更快
In [5]:
def rwh_primes2(n):
# http://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]]
In [6]:
%%timeit
rwh_primes2(10**7)
更快,使用 numpy
In [7]:
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 xrange(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)]
In [8]:
%%timeit
primesfrom2to(10**7)
In [ ]: