埃拉托斯提尼的优化筛

Question

我编写了以下代码，用于查找小于或等于n的素数。输入n = 10^7时，我的程序会崩溃。所以我想知道我的代码是否可以进一步优化。

def sieve(n):
    array = []
    for i in range(1, n):
        array.append(True)
    for i in range(0, n-1):
        for j in range(2*i+4, n+1, i+2):
            if array[j-2] != False:
                if j % (i+2) == 0:
                    array[j-2] = False
    final = []
    for i in range(0, n-1):
        if array[i] == True:
            final.append(i+2)
    return final

Answer 1

首先，我要清理内部循环中的指数。我认为+2和-2乍一看有点混乱。

def sieve(n):
    array = []
    for i in range(0, n + 1):
        array.append(True)
    for i in range(2, n + 1):
        for j in range(2*i, n+1, i):
            if array[j] != False:
                if j % (i) == 0:
                    array[j] = False
    final = []
    for i in range(0, n+1):
        if array[i] == True:
            final.append(i)
    return final

这样，数组索引直接对应于它们的数目，不需要转换。这浪费了0和1的数组空间，但我认为它更容易理解。

可能的主要优化调整循环的范围。每个非素数<= n都有一个除数<= sqrt(n)，因此我们可以将外部循环限制为range(2, math.sqrt(100) + 1)。

另外，对于非素数(因为每个非素数都有素数除数)，我们实际上不需要运行内环，我们可以添加一个if array[i] == True:来进一步减少内环的数量。

内环的范围也可以缩小。它实际上可以从i^2开始，而不是从2*i开始，因为前面的论点也适用于这里。所有小于i^2的非素数都必须有一个除数< i，因此在外部循环的早期迭代中已经设置为false。

如果应用这些更改，将得到以下代码：

def faster_sieve(n):
    array = []
    for i in range(0, n + 1):
        array.append(True)

    for i in range(2, int(math.sqrt(n)) + 1):
        if array[i] == True:
            for j in range(i*i, n + 1, i):
                if array[j] != False:
                    if j % i == 0:
                        array[j] = False
    final = []
    for i in range(2, n + 1):
        if array[i] == True:
            final.append(i)
    return final

作为比较，我们可以运行

faster_sieve(int(math.pow(10, 7)))

和

sieve(int(math.pow(10, 7)))

在我的机器上：

~ time python faster_sieve.py
python faster_sieve.py  5.44s user 0.04s system 99% cpu 5.477 total

~ time python sieve.py 
python sieve.py  32.97s user 0.03s system 99% cpu 33.003 total

更快了！

现在我们可以对python进行一些微优化，但我对python一无所知。第一步可能是将append-loop替换为更快的东西，比如array = [True] * (n + 1)，这在我的机器上节省了另一秒钟。

~ time python faster_sieve.py
python faster_sieve.py  4.46s user 0.02s system 99% cpu 4.475 total

因此，的确，您的代码可以进一步优化。

/e，我可能会补充说，在这个站点上，已经有关于筛子的python代码的良好评论了。例如，这应用了大量的python优化。

/e2

再次查看代码并查看维基百科，实际上不需要内部循环中的检查，因为j是每个循环定义的一个i的倍数，我们可以将其简化为

for j in range(i*i, n + 1, i):
    array[j] = False

它进一步优化了程序：

~ time python faster_sieve.py
python faster_sieve.py  2.61s user 0.01s system 99% cpu 2.617 total

/e3在阅读了曾傑瑞的响应并在他的实现中发现了两个bug之后，我认为还有另一个简单的优化可以添加到他的方法中。因为i在其内部循环中总是很奇怪，i * i也是如此。因此，我们可以增加内循环从i到i * 2的增量步长(i * i + i为偶数)。其结果如下：

Simon: 3.5496606826782227
Jerry (broken): 1.7031567096710205
Josay: 2.2623558044433594
Krypton: 5.4344189167022705
Krypton2: 2.5819575786590576
Jerry2: 1.396036148071289

Krypton2和Jerry2是：

def krypton2(n):
    array = [True] * (n + 1)

    for i in range(2, int(math.sqrt(n)) + 1):
        if array[i] == True:
            for j in range(i*i, n + 1, i):
                array[j] = False
    final = []
    for i in range(2, n + 1):
        if array[i] == True:
            final.append(i)
    return final


def jerry2(n):
    array = [True] * n

    result = []

    result.append(2)
    for i in range(4, n, 2):
        array[i] = False;

    for i in range(3, int(math.sqrt(n) + 1), 2):
        if array[i]:
            for j in range (i*i, n, i * 2):
                array[j] = False;

    for i in range(3, n, 2):
        if array[i]:
            result.append(i)
    return result

Answer 2

列表的

初始化

array = []
for i in range(1, n):
    array.append(True)

多次给append打电话。也许还有更好的办法。您可以使用这样的列表理解：

array = [True for _ in range(1, n)]

但是您也可以使用*操作符：

array = [True] * (n-1)

布尔检查

您可以重写：if array[j-2] != False为if array[j-2]，if array[i] == True为if array。

变量名

对于数组来说，array是一个非常糟糕的名称。给出一个能告诉这个数组内容的有趣之处的名字可能是个好主意。当您使用它来知道一个值是否是素数时，我会将它命名为is_prime，但是您可能会找到一个更好的名称。

列表理解(再次)

您可以使用列表理解重写函数的最后一部分。

    return [i + 2 for i in range(0, n-1) if is_prime[i]]

循环类似于本机

当您使用索引在Python中的数组上迭代时，通常有更好的方法来完成您想要实现的任务。我强烈建议奈德·巴奇尔德的演讲题为:像本地人一样的循环对此了解更多。

在你的例子中，你可以写：

return [i + 2 for i, p in enumerate(is_prime) if p]

更好的算法

筛的全部要点是不使用除法(或模运算符)。

我会用这样的方法：

导入数学

def sieve(n):
    is_prime = [True] * n
    is_prime[0] = is_prime[1] = False
    for i in range(2, int(math.sqrt(n)) + 1):
        if is_prime[i]:
            for j in range(i * i, n, i):
                is_prime[j] = False
    return [i for i, p in enumerate(is_prime) if p]

Answer 3

您已经收到了一些关于如何优化代码的建议，但我不能说它们中的任何一个都会让我觉得代码非常干净或简单。我并没有花很多时间试图用概要来量化对执行速度的精确影响，但在我看来，这样做是非常简单的：

def sieve(n):
    array = [True] * n

    result = []

    result.append(2)
    for i in range(4, n, 2):
        array[i] = False;

    for i in range(3, int(math.sqrt(n)), 2):
        if (array[i]):
            result.append(i)
            for j in range (i*i, n, i):
                array[j] = False;

    for i in range(int(math.sqrt(n)), n):
      if (array[i]):
        result.append(i)

这个特例是2，因为它是一个特例:唯一的偶数素数(非常特殊，许多关于素数及其性质的论文只通过排除提及2(例如，“设N是一个奇数素数”)。因此，我们首先在素数列表中添加2，然后删除2的所有倍数(如果我们想变得棘手，我们甚至不会在array中为2的倍数分配空间)。

在此之后，我们根本不需要考虑偶数，所以我们从3开始，只生成奇数作为候选。

我们只需要考虑我们极限的平方根上的数字，因为任何大于平方根的因子都必须与另一个小于平方根的因子相匹配。

对于我们找到的每一个素数，我们只需要从它平方开始就把数字划掉，因为所有比它已经被划掉的数的倍数(例如，2N已经被划出为2的倍数，3N作为3的倍数，等等)。

为了避免在找到素数后重新遍历结果列表，我们只需在找到它们时将每个元素添加到输出中。不确定这是否真的更快，但它避免了关于如何进行列表遍历的争论。

快速基准测试表明，这至少略快于我目前看到的其他任何一个：

import time
import math

def jerry(n):
    array = [True] * n

    result = []

    result.append(2)
    for i in range(4, n, 2):
        array[i] = False;

    for i in range(3, int(math.sqrt(n)), 2):
        if (array[i]):
            result.append(i)
            for j in range (i*i, n, i):
                array[j] = False;

    for i in range(int(math.sqrt(n)), n):
      if (array[i]):
        result.append(i)


def sieve(li, x):
    for i in range(2*x, len(li), x):
        li[i] = False
    return li

def simon(limit):
    li = [False, True]*(limit//2)
    li[1] = False
    li[2] = True
    for i in range(3, limit, 2):
        if li[i]:
            li = sieve(li, i)
    primes = [x for x in range(len(li)) if li[x]]
    return primes

def josay(n):
    is_prime = [True] * n
    is_prime[0] = is_prime[1] = False
    for i in range(2, int(math.sqrt(n)) + 1):
        if is_prime[i]:
            for j in range(i * i, n, i):
                is_prime[j] = False
    return [i for i, p in enumerate(is_prime) if p]

def krypton(n):
    array = []
    for i in range(0, n + 1):
        array.append(True)

    for i in range(2, int(math.sqrt(n)) + 1):
        if array[i] == True:
            for j in range(i*i, n + 1, i):
                array[j] = False
    final = []
    for i in range(2, n + 1):
        if array[i] == True:
            final.append(i)
    return final

if __name__ == '__main__':
    n = 10000000

    tests = [simon, jerry, josay, krypton]
    names = ["Simon", "Jerry", "Josay", "Krypton"]

    for i in range(0, len(tests)):
        start = time.time()
        result = tests[i](n)
        end = time.time()
        print '{}: {}'. format(names[i], (end-start))

结果：

Simon: 3.54199981689
Jerry: 2.21599984169
Josay: 3.52499985695
Krypton: 5.94400000572

注意:在发布了这些之后，@氪0n又增加了几个看起来不错的可能性，但(至少到目前为止)还没有包括在这里。