Eratosthenes筛网-Python

def sieve(upperlimit):
    # mark off all multiples of 2 so we can use 2*p as the step for the inner loop
    l = [2] + [x if x % 2 != 0 else 0 for x in range(3, upperlimit + 1)]

    for p in l:
        if p ** 2 > upperlimit:
            break
        elif p:
            for i in range(p * p, upperlimit + 1, 2 * p):
                l[i - 2] = 0
    # filter out non primes from the list, not really that important i could work with a list full of zeros as well
    return [x for x in l if x]

#1 楼

这里是计算性能改进的起点，这归因于以下各种修订：筛选\ $ 10 ^ 8 \ $以下的质数需要多长时间？

 >>> from timeit import timeit
>>> test = lambda f: timeit(lambda:f(10**8), number=1)
>>> t1 = test(sieve)
 

确切的数字将取决于计算机的速度，因此我将计算性能比率，但记录在这里：

 >>> t1
78.9875438772142
 

您的列表l的初始化花费了一半以上的时间，因此让我们尝试一种更便宜的方法。我们还要给这个数组起一个更好的名称，并在使用时将其设为布尔型数组。

def sieve2(n):
    """Return a list of the primes below n."""
    prime = [True] * n
    for p in range(3, n, 2):
        if p ** 2 > n:
            break
        if prime[p]:
            for i in range(p * p, n, 2 * p):
                prime[i] = False
    return [2] + [p for p in range(3, n, 2) if prime[p]]

优化此类函数时，始终值得保留围绕未优化版本检查优化版本的正确性：

 >>> sieve(10**6) == sieve2(10**6)
True
 

此已经运行了不到三分之一的时间：

 >>> test(sieve2) / t1
0.30390444573149544
 

我们可以通过计算更严格的循环限制来避免对p ** 2 > n进行测试。请注意，我在这里使用了n ** .5，因为它比math.sqrt(n)快一点。

def sieve3(n):
    """Return a list of the primes below n."""
    prime = [False, False, True] + [True, False] * (n // 2)
    for p in range(3, int(n ** .5) + 1, 2):
        if prime[p]:
            for i in range(p * p, n, 2 * p):
                prime[i] = False
    return [p for p in range(2, n) if prime[p]]

这与总体运行时间没有什么区别：

 >>> test(sieve3) / t1
0.2971086436068156
 

我们可以随时累积结果，而不必在最后进行单独的迭代。请注意，我已将result.append缓存在局部变量中，以避免每次在循环中查找它。

def sieve4(n):
    """Return a list of the primes below n."""
    prime = [False, False, True] + [True, False] * (n // 2)
    result = [2]
    append = result.append
    sqrt_n = (int(n ** .5) + 1) | 1    # ensure it's odd
    for p in range(3, sqrt_n, 2):
        if prime[p]:
            append(p)
            for i in range(p * p, n, 2 * p):
                prime[i] = False
    for p in range(sqrt_n, n, 2):
        if prime[p]:
            append(p)
    return result

同样，这几乎没有什么区别：

 >>> test(sieve4) / t1
0.286016401170129
 

将筛子条目设置为False时，可以使用Python的切片分配而不是循环。这看起来很浪费，因为我们先创建了一个大列表然后将其丢弃，但这避免了昂贵的for循环和相关的Python解释器开销。但有明显的改进：

 >>> test(sieve5) / t1
0.2617646381557855
 

对于数值性能有很大的改进代码，我们可以使用NumPy。

def sieve5(n):
    """Return a list of the primes below n."""
    prime = [True] * n
    result = [2]
    append = result.append
    sqrt_n = (int(n ** .5) + 1) | 1    # ensure it's odd
    for p in range(3, sqrt_n, 2):
        if prime[p]:
            append(p)
            prime[p*p::2*p] = [False] * ((n - p*p - 1) // (2*p) + 1)
    for p in range(sqrt_n, n, 2):
        if prime[p]:
            append(p)
    return result

这是sieve5的6倍以上，是原始代码的25倍以上：

 >>> test(sieve6) / t1
0.03726392181902129
 

我们可以避免为偶数分配空间，从而改善了内存局部性：

import numpy

def sieve6(n):
    """Return an array of the primes below n."""
    prime = numpy.ones(n, dtype=numpy.bool)
    prime[:2] = False
    prime[4::2] = False
    sqrt_n = int(n ** .5) + 1
    for p in range(3, sqrt_n, 2):
        if prime[p]:
            prime[p*p::2*p] = False
    return prime.nonzero()[0]

 >>> test(sieve7) / t1
0.029220096670965198
 

最后，一个实现会分别筛选\ $ 6i−1 \ $和\ $ 6i + 1 + \ $质数的实现，这要归功于Robert William Hanks：

def sieve7(n):
    """Return an array of the primes below n."""
    prime = numpy.ones(n // 2, dtype=numpy.bool)
    sqrt_n = int(n ** .5) + 1
    for p in range(3, sqrt_n, 2):
        if prime[p // 2]:
            prime[p*p // 2::p] = False
    result = 2 * prime.nonzero()[0] + 1
    result[0] = 2
    return result

这是原始sieve的约40倍：

 >>> test(sieve8) / t1
0.023447068662434022
 

#2 楼

您将所有偶数都设置为0的第一个方法不是很有效，筛分的全部目的是避免那些昂贵的模运算。请尝试以下操作：

l = range(2, upperlimit+1) # use list(range(...)) in Python 3
l[2::2] = [0] * ((len(l) - 3) // 2 + 1)

对于其他质数的筛子零点的设置，您可以执行类似的操作，但是要弄清楚有多少个零点会很复杂。在右边添加。

#3 楼

您已经正确实现了优化p ** 2 > upperlimit，但是，有一种更有效的方法可以执行此优化。该产品不是很昂贵，但也完全没有必要。...

您可以只计算一次上限的平方根，然后重用该计算值作为对p ** 2 > upperlimit的直接比较。 。考虑：

rootlimit = math.sqrt(upperlimit)
for p in l:
    if p > rootlimit:
        break;
    ....

另外，这是一件小事，但是如果您在循环中的条件循环中包含p，p或break，则无需使用continue ...只是一件小事，但是您的代码可能是：

rootlimit = math.sqrt(upperlimit)
for p in l:
    if p > rootlimit:
        break;
    if p:
        ....

#4 楼

for p in l:
    if p ** 2 > upperlimit:
        break
    elif p:
        ...

在此循环中，p通常为零，而p ** 2 > upperlimit始终为false，直到break为止。这意味着您总是要评估两个条件。如果将if p:放在第一位，则只要p不是素数，就只会评估一个条件。