使用Sundaram的Sieve，我想我打破了pure-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)

Comparasion:

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

如果你可以控制N，列出所有质数的最快方法就是预先计算它们。认真对待。预计算是一种被忽视的优化方法。

下面是一个使用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]

我在这里找到了一个纯Python 2素数生成器，在Willy Good的评论中，它比rwh2_primes快。

def primes235(limit):
yield 2; yield 3; yield 5
if limit < 7: return
modPrms = [7,11,13,17,19,23,29,31]
gaps = [4,2,4,2,4,6,2,6,4,2,4,2,4,6,2,6] # 2 loops for overflow
ndxs = [0,0,0,0,1,1,2,2,2,2,3,3,4,4,4,4,5,5,5,5,5,5,6,6,7,7,7,7,7,7]
lmtbf = (limit + 23) // 30 * 8 - 1 # integral number of wheels rounded up
lmtsqrt = (int(limit ** 0.5) - 7)
lmtsqrt = lmtsqrt // 30 * 8 + ndxs[lmtsqrt % 30] # round down on the wheel
buf = [True] * (lmtbf + 1)
for i in xrange(lmtsqrt + 1):
    if buf[i]:
        ci = i & 7; p = 30 * (i >> 3) + modPrms[ci]
        s = p * p - 7; p8 = p << 3
        for j in range(8):
            c = s // 30 * 8 + ndxs[s % 30]
            buf[c::p8] = [False] * ((lmtbf - c) // p8 + 1)
            s += p * gaps[ci]; ci += 1
for i in xrange(lmtbf - 6 + (ndxs[(limit - 7) % 30])): # adjust for extras
    if buf[i]: yield (30 * (i >> 3) + modPrms[i & 7])

我的结果:

$ time ./prime_rwh2.py 1e8
5761455 primes found < 1e8

real    0m3.201s
user    0m2.609s
sys     0m0.578s
$ time ./prime_wheel.py 1e8
5761455 primes found < 1e8

real    0m2.710s
user    0m2.469s
sys     0m0.219s

.．.在我最近的中档笔记本电脑(i5 8265U 1.6GHz)上运行Ubuntu Win 10。

这是一个mod 30轮筛，跳过倍数2,3和5。对我来说，它在2.5e9左右的时候工作得很好，那时我的笔记本电脑开始用完8G内存，需要大量交换。

我喜欢对30取余，因为它只有8个余数不是2 3 5的倍数。这允许使用移位和“&”进行乘法，除法和mod，并应该允许将一个mod 30轮的结果打包到一个字节中。我把威利的代码变成了一个分段的mod 30轮筛，以消除大N的抖动，并张贴在这里。

还有一个更快的Javascript版本，它是分段的，并使用了一个mod 210轮(没有2,3,5或7的倍数)@GordonBGood与一个深入的解释，这对我很有用。

假设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

我很惊讶居然没人提到numba。

该版本在2.47 ms±36.5µs内达到1M标记。

几年前，维基百科页面上出现了一个阿特金筛子的伪代码。这已经不存在了，参考阿特金筛似乎是一个不同的算法。一个2007/03/01版本的维基百科页面(Primer number as 2007-03-01)显示了我用作参考的伪代码。

import numpy as np
from numba import njit

@njit
def nb_primes(n):
    # Generates prime numbers 2 <= p <= n
    # Atkin's sieve -- see https://en.wikipedia.org/w/index.php?title=Prime_number&oldid=111775466
    sqrt_n = int(np.sqrt(n)) + 1

    # initialize the sieve
    s = np.full(n + 1, -1, dtype=np.int8)
    s[2] = 1
    s[3] = 1

    # put in candidate primes:
    # integers which have an odd number of
    # representations by certain quadratic forms
    for x in range(1, sqrt_n):
        x2 = x * x
        for y in range(1, sqrt_n):
            y2 = y * y
            k = 4 * x2 + y2
            if k <= n and (k % 12 == 1 or k % 12 == 5): s[k] *= -1
            k = 3 * x2 + y2
            if k <= n and (k % 12 == 7): s[k] *= -1
            k = 3 * x2 - y2
            if k <= n and x > y and k % 12 == 11: s[k] *= -1

    # eliminate composites by sieving
    for k in range(5, sqrt_n):
        if s[k]:
            k2 = k*k
            # k is prime, omit multiples of its square; this is sufficient because
            # composites which managed to get on the list cannot be square-free
            for i in range(1, n // k2 + 1):
                j = i * k2 # j ∈ {k², 2k², 3k², ..., n}
                s[j] = -1
    return np.nonzero(s>0)[0]

# initial run for "compilation" 
nb_primes(10)

时机

In[10]:
%timeit nb_primes(1_000_000)

Out[10]:
2.47 ms ± 36.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In[11]:
%timeit nb_primes(10_000_000)

Out[11]:
33.4 ms ± 373 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In[12]:
%timeit nb_primes(100_000_000)

Out[12]:
828 ms ± 5.64 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)