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]


0 False
1 False
2 True
3 True
4 False
5 True
6 False
7 True
8 False
9 False
10 False
11 True
12 False
13 True
14 False
15 False
16 False
17 True
18 False
19 True
20 False
21 False
22 False
23 True
24 False
25 False
26 False
27 False
28 False
29 True
30 False
31 True
32 False
33 False
34 False
35 False
36 False
37 True
38 False
39 False
40 False
41 True
42 False
43 True
44 False
45 False
46 False
47 True
48 False
49 False
50 False
51 False
52 False
53 True
54 False
55 False
56 False
57 False
58 False
59 True
60 False
61 True
62 False
63 False
64 False
65 False
66 False
67 True
68 False
69 False
70 False
71 True
72 False
73 True
74 False
75 False
76 False
77 False
78 False
79 True
80 False
81 False
82 False
83 True
84 False
85 False
86 False
87 False
88 False
89 True
90 False
91 False
92 False
93 False
94 False
95 False
96 False
97 True
98 False
99 False

測試一下速度


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


1 loop, best of 3: 9.77 s per loop

改進一些


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)


1 loop, best of 3: 1.06 s per loop

更快


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)


1 loop, best of 3: 766 ms per loop

更快,使用 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)


10 loops, best of 3: 78.1 ms per loop

In [ ]: