如何加速 Python 中的埃拉托色尼筛法列表生成器

Question

我的问题直接来源于CS circles网站。它是此页面底部的最后一个问题，名为“起飞准备”。基本情况是，他们想要一个长度为1,000,001的列表，其中每个项目的索引如果是质数，则该项目的索引为True，否则为False。

例如，isPrime [13]为True。 isPrime [14]为False。

列表'isPrime'的前几个小部分如下：

isPrime = [False, False, True, True, False, True, False, True, False, False, ...]

我的问题是他们规定了7秒的时间限制。我有一个下面的工作代码，为了调试目的，它使用更低的数字101。当我将它增加到他们要求的1000001列表长度时，它需要太长时间，几分钟后我不得不结束程序。

n = 101
c = 2
isPrime = [False,False]
for i in range(2,n):
    isPrime.append(i)

def ifInt(isPrime):
    for item in isPrime:
        if type(item) == int:
            return item
for d in range(2,n):
    if c != None:
        for i in range(c,n,c):
            isPrime[i] = False
        isPrime[c] = True
        c = ifInt(isPrime)

这里有一个我找到的链接，它的运行时间更快，但只输出素数列表，而不是长度为n的列表，其中list[n]返回素数为True，否则为False。

我不确定这段代码是否真的能在不到7秒的时间内创建一个长度为1,000,001的素数列表，但它是我研究中找到的最相关的东西。

def primes_sieve1(limit):
    limitn = limit+1
    primes = dict()
    for i in range(2, limitn): primes[i] = True

    for i in primes:
        factors = range(i,limitn, i)
        for f in factors[1:]:
            primes[f] = False
    return [i for i in primes if primes[i]==True]

CS圆圈似乎很常用，但我找不到Python的可用版本。希望有人能轻松解决这个问题。

这个问题与这个问题不同，因为我不仅想快速创建质数列表，还想创建从0到n的所有正整数列表，其中被标记为True的是素数，被标记为False的是非素数。

- adhdj

---

以下是需要翻译的内容：https://dev59.com/snI95IYBdhLWcg3w-DH0答案：这个问题有很多算法可以解决。其中一种常见的算法是埃拉托色尼筛法（Sieve of Eratosthenes）。它的基本思想是从小到大枚举素数并将其的倍数标记为非素数。这样，最终没有被标记的数就是素数。以下是 Python 的实现代码： - Untitled123

请注意，列出质数与您所描述的True / False数组完全等效，并且可能任何列出质数的程序只需要进行非常小的更改即可获得布尔数组。 - Untitled123

@adhdj，你所链接的答案中的代码，包括打印输出，比如print list(sieve.primes_sieve2(10000000)) 在不到7秒的时间内运行。它比你的问题大10倍，并在内部生成了问题中所需的列表。你所要做的就是将其取出即可。 - pvg

好的，我会试着调整一下，但我的问题不是在7秒内获取质数列表，而是在7秒内将该列表按照上述格式排列。不过我认为我已经走在了正确的道路上。@pvg - adhdj

@adhdj 我不确定如何以更清晰的方式重复这一点。你链接的答案_生成了所讨论的列表_。 - pvg

Answer 1

我意识到在SO上有很多针对素数筛法的优化，但其他人很少解释它们，因此对初学者或第一次创建该算法的人来说，这些优化很难接近。这里所有的解决方案都是用Python编写的，以达到相同的速度和最优化。这些解决方案将逐渐变得更快，更复杂。 :)

基本解决方案

def primesVanilla(n):
    r = [True] * n
    r[0] = r[1] = False
    for i in xrange(n):
        if r[i]:
            for j in xrange(i+i, n, i):
                r[j] = False
    return r

这是一个非常直观的筛法实现。在继续之前，请确保您理解上面发生的事情。唯一需要注意的是，您从i+i（而不是从i）开始标记非质数，但这是相当明显的。（因为您假设i本身是一个质数）。为了使测试公平，所有数字将针对列表最多2500万个进行。

real    0m7.663s  
user    0m7.624s  
sys     0m0.036s  

小改进1（平方根）：

我会尝试按照直截了当到不那么直截了当的顺序进行排序。请注意，我们不需要迭代到n，而是只需要到n的平方根。原因是在n以下的任何合数，必须具有小于或等于n的平方根的质因子。当您手动筛选时，您会注意到所有“未筛选”的大于n的平方根的数字默认为质数。

另一个要注意的是，当平方根是整数时，必须要小心一点，所以在这种情况下应该加上一，以覆盖它。例如，在n = 49时，您希望循环到7（包括7），否则您可能会得出49是质数的结论。

def primes1(n):
    r = [True] * n
    r[0] = r[1] = False
    for i in xrange(int(n**0.5+1)):
        if r[i]:
            for j in xrange(i+i, n, i):
                r[j] = False
    return r

real    0m4.615s
user    0m4.572s
sys     0m0.040s

请注意，这样做要快得多。当您考虑到它时，您只需要循环到平方根，所以现在需要2500万个顶级迭代的东西现在只需要5000个顶级迭代。

小改进2（在内部循环中跳过）：

观察一下，在内部循环中，我们可以从i * i开始而不是从i + i开始。这与平方根的情况非常相似，但大的想法是i和i * i之间的任何复合数已经被较小的质数标记过了。

def primes2(n):
    r = [True] * n
    r[0] = r[1] = False
    for i in xrange(int(n**0.5+1)):
        if r[i]:
            for j in xrange(i*i, n, i):
                r[j]=False
    return r

real    0m4.559s
user    0m4.500s
sys     0m0.056s

嗯，有点令人失望。但是，嘿，它仍然更快了。

较大的改进3（甚至跳过）：

这里的想法是我们可以预先标记所有偶数索引，然后在主循环中每2次迭代跳过一次。之后，我们可以从3开始外部循环，内部循环可以跳过2 * i。 (因为按i进行循环意味着它将是偶数，(i+i) (i+i+i+i)等等)

def primes3(n):
    r = [True] * n
    r[0] = r[1] = False
    for i in xrange(4,n,2):
        r[i] = False    
    for i in xrange(3, int(n**0.5+1), 2):
        if r[i]:
            for j in xrange(i*i, n, 2*i):
                r[j] = False
    return r

real    0m2.916s
user    0m2.872s
sys     0m0.040s

Cool Improvements 4（wim的想法）：

这个解决方案是一个相当高级的技巧。切片赋值比循环更快，因此它使用Python的切片符号：r[begin:end:skip]

def primes4(n):
    r = [True] * n
    r[0] = r[1] = False 
    r[4::2] = [False] * len(r[4::2])
    for i in xrange(3, int(1 + n**0.5), 2):
        if r[i]:
            r[i*i::2*i] = [False] * len(r[i*i::2*i])
    return r

10 loops, best of 3: 1.1 sec per loop

轻微改进 5

请注意，Python在计算长度时会重新切片r[4::2]，因此这会花费相当多的额外时间，因为我们所需的只是计算长度。不过，我们确实使用了一些不好的数学方法来实现这一点。

def primes5(n):
    r = [True] * n
    r[0] = r[1] = False 
    r[4::2] = [False] * ((n+1)/2-2)
    for i in xrange(3, int(1 + n**0.5), 2):
        if r[i]:
            r[i*i::2*i] = [False] * ((n+2*i-1-i*i)/(2*i))
    return r

10 loops, best of 3: 767 msec per loop

作业加速（Padraic Cunningham）：

注意，我们分配一个包含全部 True 值的数组，然后将一半（偶数位）设置为 False。实际上，我们可以直接使用交替的布尔数组。

def primes6(n):
    r = [False, True] * (n//2) + [True]
    r[1], r[2] = False, True
    for i in xrange(3, int(1 + n**0.5), 2):
        if r[i]:
            r[i*i::2*i] = [False] * ((n+2*i-1-i*i)/(2*i))
    return r

10 loops, best of 3: 717 msec per loop

不要引用我的话，但我认为在没有一些复杂的数学方法的情况下，对于这个最新版本没有明显的改进。我尝试过一个聪明的方法，但并没有变得更快，即注意到除了2和3之外的素数必须是6k + 1或6k-1的形式。（请注意，如果它是6k，则可被6整除，6k + 2 | 2，6k + 3 | 3，6k + 4 | 2，6k + 5与模6同余于-1。这表明我们可以每次跳过6并检查两侧。但无论是因为我实现得不好还是Python内部的原因，我都无法找到任何有意义的速度提升。 :(

- Untitled123

---

+1 很好。我一直认为Python中的列表切片应该是视图而不是副本。我很遗憾这在Python3中没有实现。 - wim

创建一个切片仅仅为了找到它的长度是非常浪费时间和内存的，当我们可以很容易地计算出正确的长度。请参见我的布尔素数筛选器。 - PM 2Ring

我仍然遇到“ValueError：尝试将大小为49999的序列分配给大小为49998的扩展切片”的错误。 - Padraic Cunningham

你可能想把那些 / 改成 //，以在 Python 3 中强制执行整数除法。 - PM 2Ring

xrange 也是 Python 2.7 的遗迹，我的大脑告诉我要转向3，但我的内心告诉我不要。 - Untitled123

Answer 2

我看到的第一件事是，您生成初始列表的方式(循环和追加)是低效且不必要的。您可以直接添加列表，而无需逐个循环和追加。

我看到的第二件事是，您正在进行的类型检查是不必要的，该函数调用是昂贵的，您可以重构以完全避免它。

最后，我认为您在任何筛子实现中都能得到的“大事情”就是利用切片赋值。您应该一次性删除所有因子，而不是循环。例如：

from math import sqrt

def primes(n):
    r = [True] * n
    r[0] = r[1] = False
    r[4::2] = [False] * len(r[4::2])
    for i in xrange(int(1 + sqrt(n))):
        if r[i]:
            r[3*i::2*i] = [False] * len(r[3*i::2*i])
    return r

注意，我还有一些其他的技巧：

• 通过立即划掉偶数来避免一半的工作。
• 只需要迭代到长度的平方根就足够了。

在我的破旧低配MacBook上，这段代码可以在大约75毫秒内生成1,000,001列表：

>>> timeit primes(1000001)
10 loops, best of 3: 75.4 ms per loop

- wim

---

我认为你应该使用range而不是enumerate，因为它会复制列表。 - simonzack

枚举不行，但切片可以。 - simonzack

啊，明白了。好主意，我会修复它的。 - wim

我认为即使交叉运行两个程序，也不会节省任何工作量。目前的措辞有些欺骗性吗？此外，已经有非常好的答案来实现疯狂快速的素数筛和优化技巧了...为什么不链接到其中一个呢？ - Untitled123

如果您能在保持简单性的同时提高效率，我会非常感兴趣看看 :) - wim

Answer 3

在Python2和3中，一些时间显示wim的方法明显更快，可以通过创建列表的方式进一步进行轻微优化：

def primes_wim_opt(n):
    r = [False, True] * (n // 2)
    r[0] = r[1] = False
    r[2] = True
    for i in xrange(int(1 + n ** .5)):
        if r[i]:
            r[3*i::2*i] = [False] * len(r[3*i::2*i])
    return r

Python2的时间记录：

In [9]: timeit primesVanilla(100000)
10 loops, best of 3: 25.7 ms per loop

In [10]: timeit primes_wim(100000)
100 loops, best of 3: 3.59 ms per loop

In [11]: timeit primes1(100000)
100 loops, best of 3: 14.8 ms per loop

In [12]: timeit primes_wim_opt(100000)
100 loops, best of 3: 2.18 ms per loop

In [13]: timeit primes2(100000)
100 loops, best of 3: 14.7 ms per loop

In [14]: primes_wim(100000) ==  primes_wim_opt(100000) ==  primes(100000) == primesVanilla(100000) == primes2(100000)
Out[14]: True

使用相同的函数，只是将其更改为range()的python3时间：

In [76]: timeit primesVanilla(100000)
10 loops, best of 3: 22.3 ms per loop

In [77]: timeit primes_wim(100000)
100 loops, best of 3: 2.92 ms per loop

In [78]: timeit primes1(100000)
100 loops, best of 3: 10.9 ms per loop

In [79]: timeit primes_wim_opt(100000)
1000 loops, best of 3: 1.88 ms per loop

In [80]: timeit primes2(100000)
100 loops, best of 3: 10.3 ms per loop
In [81]: primes_wim(100000) ==  primes_wim_opt(100000) ==  primes(100000) == primesVanilla(100000) == primes2(100000)
Out[80]: True

可以通过使用range/xrange的长度而不是切片来进一步优化：

def primes_wim_opt(n):
    is_odd = n % 2 & 1    
    r = [False, True] * (n // 2 + is_odd)
    r[0] = r[1] = False
    r[2] = True
    for i in range(int(1 + n ** .5)):
        if r[i]:
            r[3*i::2*i] = [False] * len(range(3*i,len(r), 2 * i))
    return r

Python3可以大幅提高IT技术的效率：

In [16]: timeit primes_wim_opt_2(100000)
1000 loops, best of 3: 1.38 ms per loop

对于使用xrange的Python2也是如此：

In [10]: timeit  primes_wim_opt_2(100000)
1000 loops, best of 3: 1.60 ms per loop

使用(((n - 3 * i) // (2 * i)) + 1)也可以起作用：

def primes_wim_opt_2(n):
    is_odd = n % 2 & 1
    r = [False, True] * ((n // 2) + is_odd)
    r[0] = r[1] = False
    r[2] = True
    for i in range(int(1 + n ** .5)):
        if r[i]:
            r[3*i::2*i] = [False] * (((n - 3 * i) // (2 * i)) + 1)
    return r

哪一个速度稍微快一点：

In [12]: timeit primes_wim_opt_2(100000)
1000 loops, best of 3: 1.32 ms per loop

In [6]: timeit primes5(100000)
100 loops, best of 3: 2.47 ms per loop

您也可以从3开始，步长为2：

def primes_wim_opt_2(n):
    r = [False, True] * (n // 2)
    r[0] = r[1] = False
    r[2] = True
    for i in range(3, int(1 + n ** .5),2):
        if r[i]:
            r[3*i::2*i] = [False] * (((n - 3 * i) // (2 * i)) + 1)
    return r

哪个更快：

In [2]: timeit primes_wim_opt_2(100000)
1000 loops, best of 3: 1.10 ms per loop

Python2：

In [2]: timeit primes_wim_opt_2(100000)
1000 loops, best of 3: 1.29 ms per loop

- Padraic Cunningham

---

你也成为了计算的受害者：）对于1000001的切片长度问题。 - Untitled123

你介意我在我的回答中加入[False,True]技巧吗？可以提高10%的速度 :) - Untitled123

使用range/xrange的长度而不是切片是个好主意。这比使用预先计算的大分数要少出错得多。 - wim