如何加速埃拉托色尼筛法 python 列表生成器
How to speed up Sieve of Eratosthenes python list generator
我的问题直接来自 CS 圈子网站。这是 this 页面底部的最后一个问题 'Primed for Takeoff'。基本的 运行down 是他们想要一个 1,000,001 长度的列表，其中如果索引是素数，则每个项目的索引为 True，而如果不是素数，则每个项目的索引为 False。
因此，例如，isPrime[13] 为真。 isPrime[14] 为假。
列表的第一部分 'isPrime' 看起来像：
isPrime = [False, False, True, True, False, True, False, True, False, False, ...]
我的问题是他们有 7 秒的时间限制。为了调试目的，我在下面有一个较小的数字 101 的工作代码。当我将其增加到他们要求的 1,000,001 列表长度时，它花费的时间太长了，几分钟后我结束了程序。
这是我的工作（长度为 101），非常慢的代码：
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)
然后就是我找到的这个 here。它有更快的 运行 时间，但只输出素数列表，而不是 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]
print primes_sieve1(101)
CS 圈子似乎很常用，但我找不到 Python 的可用版本。希望这对某人来说是一个简单的解决方案。
这个问题与this one不同，因为我不只是试图快速创建一个素数列表，而是创建一个从 0 到 n 的所有正整数的列表，这些正整数被 True 和非素数标记为素数假的。
我首先看到的是生成初始列表（循环和附加）的方式效率低下且不必要。您可以只 添加 列表而不是循环和附加 per-element。
我看到的第二件事是你正在做的 type-checking 是不必要的，那个函数调用很昂贵，你可以重构来完全避免这种情况。
最后，我认为 "big thing" 您可以在任何筛子实现中获得的是利用切片分配。您应该一次性划掉所有因素，而不是循环。示例：
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
注意我还有一些其他技巧：
- 立即划掉偶数可以避免一半的工作。
- 只需要迭代到长度的 sqrt
在我性能低下的 macbook 上，这段代码可以在大约 75 毫秒内生成 1,000,001 个列表：
>>> timeit primes(1000001)
10 loops, best of 3: 75.4 ms per loop
我意识到 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
这是 Sieve 的一个非常简单的实现。在继续之前，请确保您了解上面发生的事情。唯一需要注意的一点是，您开始在 i+i 而不是 i 处标记 not-primes，但这是相当明显的。 （因为您假设 i 本身是素数）。为了使测试公平，所有数字将用于列表最多 2500 万.
real 0m7.663s
user 0m7.624s
sys 0m0.036s
次要改进 1（平方根）：
我将尝试根据 straight-forward 对它们进行排序以减少 straight-forward 更改。观察到我们不需要迭代到 n，而只需要达到 n 的平方根。原因是任何 n 以下的合数都必须有一个小于或等于 n 的平方根的质因数。当您手动筛选时，您会注意到 n 的平方根上的所有“未筛选”数字默认都是素数。
另一个注意事项是当平方根变成整数时你必须小心一点，所以在这种情况下你应该加一个以覆盖它。 IE，在 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
很酷的改进 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
作业speed-up（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 整除为 -1 mod 6。这表明我们可以跳过每次增加 6 并检查两侧。无论是我这边的实施不当，还是 python 内部结构，我都找不到任何有意义的速度提升。:(
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
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 的 len 而不是切片来进一步优化：
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 它敲掉了一大块：
In [16]: timeit primes_wim_opt_2(100000)
1000 loops, best of 3: 1.38 ms per loop
python2 使用 xrange:
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
我的问题直接来自 CS 圈子网站。这是 this 页面底部的最后一个问题 'Primed for Takeoff'。基本的 运行down 是他们想要一个 1,000,001 长度的列表，其中如果索引是素数，则每个项目的索引为 True，而如果不是素数，则每个项目的索引为 False。
因此，例如，isPrime[13] 为真。 isPrime[14] 为假。
列表的第一部分 'isPrime' 看起来像：
isPrime = [False, False, True, True, False, True, False, True, False, False, ...]
我的问题是他们有 7 秒的时间限制。为了调试目的，我在下面有一个较小的数字 101 的工作代码。当我将其增加到他们要求的 1,000,001 列表长度时，它花费的时间太长了，几分钟后我结束了程序。 这是我的工作（长度为 101），非常慢的代码：
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)
然后就是我找到的这个 here。它有更快的 运行 时间，但只输出素数列表，而不是 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]
print primes_sieve1(101)
CS 圈子似乎很常用，但我找不到 Python 的可用版本。希望这对某人来说是一个简单的解决方案。
这个问题与this one不同，因为我不只是试图快速创建一个素数列表，而是创建一个从 0 到 n 的所有正整数的列表，这些正整数被 True 和非素数标记为素数假的。
我首先看到的是生成初始列表（循环和附加）的方式效率低下且不必要。您可以只 添加 列表而不是循环和附加 per-element。
我看到的第二件事是你正在做的 type-checking 是不必要的，那个函数调用很昂贵，你可以重构来完全避免这种情况。
最后，我认为 "big thing" 您可以在任何筛子实现中获得的是利用切片分配。您应该一次性划掉所有因素，而不是循环。示例：
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
注意我还有一些其他技巧：
- 立即划掉偶数可以避免一半的工作。
- 只需要迭代到长度的 sqrt
在我性能低下的 macbook 上，这段代码可以在大约 75 毫秒内生成 1,000,001 个列表：
>>> timeit primes(1000001)
10 loops, best of 3: 75.4 ms per loop
我意识到 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
这是 Sieve 的一个非常简单的实现。在继续之前，请确保您了解上面发生的事情。唯一需要注意的一点是，您开始在 i+i 而不是 i 处标记 not-primes，但这是相当明显的。 （因为您假设 i 本身是素数）。为了使测试公平，所有数字将用于列表最多 2500 万.
real 0m7.663s
user 0m7.624s
sys 0m0.036s
次要改进 1（平方根）：
我将尝试根据 straight-forward 对它们进行排序以减少 straight-forward 更改。观察到我们不需要迭代到 n，而只需要达到 n 的平方根。原因是任何 n 以下的合数都必须有一个小于或等于 n 的平方根的质因数。当您手动筛选时，您会注意到 n 的平方根上的所有“未筛选”数字默认都是素数。
另一个注意事项是当平方根变成整数时你必须小心一点，所以在这种情况下你应该加一个以覆盖它。 IE，在 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
很酷的改进 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
作业speed-up（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 整除为 -1 mod 6。这表明我们可以跳过每次增加 6 并检查两侧。无论是我这边的实施不当，还是 python 内部结构，我都找不到任何有意义的速度提升。:(
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
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 的 len 而不是切片来进一步优化：
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 它敲掉了一大块：
In [16]: timeit primes_wim_opt_2(100000)
1000 loops, best of 3: 1.38 ms per loop
python2 使用 xrange:
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