这是我能想到的最好的算法。

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

当前回答

我在这里找到了一个纯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与一个深入的解释,这对我很有用。

其他回答

到目前为止,我尝试过的最快的方法是基于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

关于加速的解释,请看下面的答案。

我猜最快的方法是在代码中硬编码质数。

因此,为什么不编写一个缓慢的脚本,生成另一个源文件,其中包含所有数字,然后在运行实际程序时导入该源文件呢?

当然,只有当你在编译时知道N的上限时,这才有效,但这是(几乎)所有项目欧拉问题的情况。

 

PS:我可能错了,虽然解析源的硬连接质数比计算它们要慢,但据我所知,Python是从编译的.pyc文件运行的,所以在这种情况下,读取一个包含所有质数到N的二进制数组应该是非常快的。

我可能迟到了,但必须为此添加自己的代码。它使用大约n/2的空间,因为我们不需要存储偶数,我还使用bitarray python模块,进一步大幅减少内存消耗,并允许计算所有高达1,000,000,000的质数

from 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

这是在64bit 2.4GHZ MAC OSX 10.8.3上运行的

我对这个问题反应迟钝,但这似乎是一个有趣的练习。我使用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]

对于足够大的N,真正最快的解决方案是下载一个预先计算的质数列表,将其存储为元组,并执行如下操作:

for pos,i in enumerate(primes):
    if i > N:
        print primes[:pos]

如果只有N >个质数[-1],则计算更多的质数并将新列表保存在代码中,以便下次同样快。

要跳出思维定势。