列出N以下所有质数的最快方法

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

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

当前回答

这个算法很快，但它有一个严重的缺陷:

>>> sorted(get_primes(530))
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73,
79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163,
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 527, 529]
>>> 17*31
527
>>> 23*23
529

您假设numbers.pop()将返回集合中最小的数字，但这根本不能保证。集合是无序的，pop()删除并返回任意元素，因此不能使用它从剩余数字中选择下一个质数。

2010-01-15 00:14:04

其他回答

在Pure Python中最快的质数筛分:

from itertools import compress

def half_sieve(n):
    """
    Returns a list of prime numbers less than `n`.
    """
    if n <= 2:
        return []
    sieve = bytearray([True]) * (n // 2)
    for i in range(3, int(n ** 0.5) + 1, 2):
        if sieve[i // 2]:
            sieve[i * i // 2::i] = bytearray((n - i * i - 1) // (2 * i) + 1)
    primes = list(compress(range(1, n, 2), sieve))
    primes[0] = 2
    return primes

我优化了埃拉托色尼筛子的速度和内存。

基准

from time import clock
import platform

def benchmark(iterations, limit):
    start = clock()
    for x in range(iterations):
        half_sieve(limit)
    end = clock() - start
    print(f'{end/iterations:.4f} seconds for primes < {limit}')

if __name__ == '__main__':
    print(platform.python_version())
    print(platform.platform())
    print(platform.processor())
    it = 10
    for pw in range(4, 9):
        benchmark(it, 10**pw)

输出

>>> 3.6.7
>>> Windows-10-10.0.17763-SP0
>>> Intel64 Family 6 Model 78 Stepping 3, GenuineIntel
>>> 0.0003 seconds for primes < 10000
>>> 0.0021 seconds for primes < 100000
>>> 0.0204 seconds for primes < 1000000
>>> 0.2389 seconds for primes < 10000000
>>> 2.6702 seconds for primes < 100000000

2019-05-24 01:51:04

对于最快的代码，numpy解决方案是最好的。不过，出于纯粹的学术原因，我发布了我的纯python版本，它比上面发布的食谱版本快不到50%。由于我将整个列表放在内存中，所以需要足够的空间来容纳所有内容，但它的可伸缩性似乎相当好。

def daniel_sieve_2(maxNumber):
    """
    Given a number, returns all numbers less than or equal to
    that number which are prime.
    """
    allNumbers = range(3, maxNumber+1, 2)
    for mIndex, number in enumerate(xrange(3, maxNumber+1, 2)):
        if allNumbers[mIndex] == 0:
            continue
        # now set all multiples to 0
        for index in xrange(mIndex+number, (maxNumber-3)/2+1, number):
            allNumbers[index] = 0
    return [2] + filter(lambda n: n!=0, allNumbers)

结果是:

>>>mine = timeit.Timer("daniel_sieve_2(1000000)",
...                    "from sieves import daniel_sieve_2")
>>>prev = timeit.Timer("get_primes_erat(1000000)",
...                    "from sieves import get_primes_erat")
>>>print "Mine: {0:0.4f} ms".format(min(mine.repeat(3, 1))*1000)
Mine: 428.9446 ms
>>>print "Previous Best {0:0.4f} ms".format(min(prev.repeat(3, 1))*1000)
Previous Best 621.3581 ms

2010-01-15 07:41:25

我知道比赛已经结束好几年了。.．.

尽管如此，这是我对纯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

2013-09-22 13:35:58

这是问题解的一种变化应该比问题本身更快。它使用埃拉托色尼的静态筛，没有其他优化。

from typing import List

def list_primes(limit: int) -> List[int]:
    primes = set(range(2, limit + 1))
    for i in range(2, limit + 1):
        if i in primes:
            primes.difference_update(set(list(range(i, limit + 1, i))[1:]))
    return sorted(primes)

>>> list_primes(100)
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

2020-08-15 16:35:56

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

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

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

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

2010-01-25 01:21:22

列出N以下所有质数的最快方法

推荐文章

最新文章

标签