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

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、3和5的倍数。尽管如此,在N<10^9时,它实际上比@Robert William Hanks的优解rwh_primes2和rwh_primes1要慢。通过使用大于1.5* 10^8的ctypes.c_ushort筛分数组,可以在某种程度上适应内存限制。

10^6

$ python -mtimeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp. primesieveseq (1000000)" 10个循环,最好的3:46.7毫秒每循环

import primeSieveSpeedComp (primeSieveSpeedComp) “primeSieveSpeedComp.rwh_primes1(1000000)”10个循环,最好的3:43.2 每回路Msec $ python -m timeit -s"import primeSieveSpeedComp" “primeSieveSpeedComp.rwh_primes2(1000000)”10圈,最好成绩是3:34.5 每回路Msec

10^7

$ python -mtimeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp. primesieveseq (10000000)" 10个循环,最好是3:530毫秒每循环

import primeSieveSpeedComp (primeSieveSpeedComp) “primeSieveSpeedComp.rwh_primes1(10000000)”10圈,3:494的最佳成绩 每回路Msec $ python -m timeit -s"import primeSieveSpeedComp" “primeSieveSpeedComp.rwh_primes2(10000000)”10圈,最好的3:375 每回路Msec

10^8

$ python -mtimeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp. primesieveseq (100000000)" 10圈,最好的3:5.55秒每圈

import primeSieveSpeedComp (primeSieveSpeedComp) “primeSieveSpeedComp.rwh_primes1(100000000)”10圈,最好成绩是3:5.33 秒/循环 $ python -m timeit -s"import primeSieveSpeedComp" “primeSieveSpeedComp.rwh_primes2(100000000)”10圈,最好的3:3.95 秒/循环

10^9

$ python -mtimeit -s"import primeSieveSpeedComp" "primeSieveSpeedComp. primesieveseq (1000000000)" 10圈,最好的3圈:每圈61.2秒

$ python -mtimeit -n 3 -s"import primeSieveSpeedComp" “primeSieveSpeedComp.rwh_primes1(1000000000)”3圈,最好的3:97.8 秒/循环 $ python -m timeit -s"import primeSieveSpeedComp" “primeSieveSpeedComp.rwh_primes2(1000000000)”10个循环,3个最好: 每循环41.9秒

您可以将下面的代码复制到ubuntu primeSieveSpeedComp中以查看此测试。

def primeSieveSeq(MAX_Int):
    if MAX_Int > 5*10**8:
        import ctypes
        int16Array = ctypes.c_ushort * (MAX_Int >> 1)
        sieve = int16Array()
        #print 'uses ctypes "unsigned short int Array"'
    else:
        sieve = (MAX_Int >> 1) * [False]
        #print 'uses python list() of long long int'
    if MAX_Int < 10**8:
        sieve[4::3] = [True]*((MAX_Int - 8)/6+1)
        sieve[12::5] = [True]*((MAX_Int - 24)/10+1)
    r = [2, 3, 5]
    n = 0
    for i in xrange(int(MAX_Int**0.5)/30+1):
        n += 3
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 2
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 1
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 2
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 1
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 2
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 3
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
        n += 1
        if not sieve[n]:
            n2 = (n << 1) + 1
            r.append(n2)
            n2q = (n2**2) >> 1
            sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
    if MAX_Int < 10**8:
        return [2, 3, 5]+[(p << 1) + 1 for p in [n for n in xrange(3, MAX_Int >> 1) if not sieve[n]]]
    n = n >> 1
    try:
        for i in xrange((MAX_Int-2*n)/30 + 1):
            n += 3
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 2
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 1
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 2
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 1
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 2
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 3
            if not sieve[n]:
                r.append((n << 1) + 1)
            n += 1
            if not sieve[n]:
                r.append((n << 1) + 1)
    except:
        pass
    return r

其他回答

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

这里有一个来自Python Cookbook的非常简洁的示例——该URL的最快版本是:

import itertools
def erat2( ):
    D = {  }
    yield 2
    for q in itertools.islice(itertools.count(3), 0, None, 2):
        p = D.pop(q, None)
        if p is None:
            D[q*q] = q
            yield q
        else:
            x = p + q
            while x in D or not (x&1):
                x += p
            D[x] = p

这就给出了

def get_primes_erat(n):
  return list(itertools.takewhile(lambda p: p<n, erat2()))

在shell提示符(正如我喜欢做的那样)中测量这段代码在pri.py中,我观察到:

$ python2.5 -mtimeit -s'import pri' 'pri.get_primes(1000000)'
10 loops, best of 3: 1.69 sec per loop
$ python2.5 -mtimeit -s'import pri' 'pri.get_primes_erat(1000000)'
10 loops, best of 3: 673 msec per loop

所以看起来食谱解决方案的速度是原来的两倍多。

如果你接受itertools,但不接受numpy,这里有一个针对Python 3的rwh_primes2的改编版本,它在我的机器上运行速度大约是原来的两倍。唯一的实质性变化是使用bytearray而不是列表来表示布尔值,并使用压缩而不是列表推导来构建最终列表。(如果可以的话,我会把这句话作为moarningsun之类的评论。)

import itertools
izip = itertools.zip_longest
chain = itertools.chain.from_iterable
compress = itertools.compress
def rwh_primes2_python3(n):
    """ Input n>=6, Returns a list of primes, 2 <= p < n """
    zero = bytearray([False])
    size = n//3 + (n % 6 == 2)
    sieve = bytearray([True]) * size
    sieve[0] = False
    for i in range(int(n**0.5)//3+1):
      if sieve[i]:
        k=3*i+1|1
        start = (k*k+4*k-2*k*(i&1))//3
        sieve[(k*k)//3::2*k]=zero*((size - (k*k)//3 - 1) // (2 * k) + 1)
        sieve[  start ::2*k]=zero*((size -   start  - 1) // (2 * k) + 1)
    ans = [2,3]
    poss = chain(izip(*[range(i, n, 6) for i in (1,5)]))
    ans.extend(compress(poss, sieve))
    return ans

比较:

>>> timeit.timeit('primes.rwh_primes2(10**6)', setup='import primes', number=1)
0.0652179726976101
>>> timeit.timeit('primes.rwh_primes2_python3(10**6)', setup='import primes', number=1)
0.03267321276325674

and

>>> timeit.timeit('primes.rwh_primes2(10**8)', setup='import primes', number=1)
6.394284538007014
>>> timeit.timeit('primes.rwh_primes2_python3(10**8)', setup='import primes', number=1)
3.833829450302801

我可能迟到了,但必须为此添加自己的代码。它使用大约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上运行的

这是你和别人比较的方式。

# You have to list primes upto n
nums = xrange(2, n)
for i in range(2, 10):
    nums = filter(lambda s: s==i or s%i, nums)
print nums

这么简单……