埃拉托色尼筛 - 寻找素数 Python答案

【问题标题】：Sieve of Eratosthenes - Finding Primes Python埃拉托色尼筛 - 寻找素数 Python
【发布时间】：2011-04-25 18:37:54
【问题描述】：

澄清一下，这不是作业问题:)

我想为我正在构建的数学应用程序找到素数，并遇到了Sieve of Eratosthenes 方法。

我已经用 Python 编写了它的实现。但这非常慢。比如说，如果我想找到所有小于 200 万的素数。它需要> 20分钟。（我在这一点上停止了它）。如何加快速度？

def primes_sieve(limit):
    limitn = limit+1
    primes = range(2, limitn)

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

print primes_sieve(2000)

更新： 我最终对这段代码进行了分析，发现从列表中删除一个元素花了很多时间。考虑到它必须遍历整个列表（最坏情况）才能找到元素然后删除它然后重新调整列表（也许还有一些副本？），这是可以理解的。无论如何，我把字典列表扔掉了。我的新实现 -

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(2000000)

【问题讨论】：

这里有一个类似的问题stackoverflow.com/questions/2897297，您可能会觉得有用。
检查that 答案。
stackoverflow.com/questions/2068372/…
@Srikar：您可以只迭代到 limit 的平方根，而不是迭代到 limit 的平方根，因为字典中的任何合数都会比 limit 的平方根少一个因子。
对range 使用step 参数非常棒。 factors 用词不当，应为 multiples。

标签： python math primes sieve-of-eratosthenes

【解决方案1】：

你没有完全实现正确的算法：

在您的第一个示例中，primes_sieve 不维护要触发/取消设置的素数标志列表（如算法中所示），而是不断调整整数列表的大小，这非常昂贵：从list 需要将所有后续项目向下移动一个。

在第二个示例中，primes_sieve1 维护一个 dictionary 素数标志，这是朝着正确方向迈出的一步，但它以未定义的顺序迭代字典，并多余地删除了因数（而不是算法中的仅质数因数）。您可以通过对键进行排序并跳过非素数来解决此问题（这已经使其速度提高了一个数量级），但直接使用列表仍然更有效。

正确的算法（使用列表而不是字典）看起来像：

def primes_sieve2(limit):
    a = [True] * limit                          # Initialize the primality list
    a[0] = a[1] = False

    for (i, isprime) in enumerate(a):
        if isprime:
            yield i
            for n in range(i*i, limit, i):     # Mark factors non-prime
                a[n] = False

（请注意，这还包括在素数的正方形 (i*i) 而不是其双精度数处开始非素数标记的算法优化。）

【讨论】：

又一个优化，你的xrange(i*i,limit,i)的步长可以设为2*i
我喜欢你对埃拉托色尼筛法的简洁实现。 :) 但是，我遇到了 OverflowError: Python int too large to convert to C long。我将 xrange(i*i, limit, i) 更改为 xrange(i, limit, i)。感谢分享此代码 sn-p！
@st0le: 不，不能设置步长2*i。刚试了一下。它产生 14 作为素数。
@Mark，很抱歉我没有真正完整地解释它。通过使用i=2 以i 的步长进行迭代来消除所有偶数，但其余的您可以使用2*i。事实上，在我的实现中，我使用了一半的布尔值，因为我不存储偶数，而是使用简单的mod 2。你可以在这里找到我的 Java 实现，它使用更少（1/8）的内存。 HERE
+1，只是一个小细节，如果你在初始化时使用[False] * 2 + [True] * (limit-2)，你可以避免IndexError on passing number

【解决方案2】：

def eratosthenes(n):
    multiples = []
    for i in range(2, n+1):
        if i not in multiples:
            print (i)
            for j in range(i*i, n+1, i):
                multiples.append(j)

eratosthenes(100)

【讨论】：

我会使用集合而不是列表，以加快成员资格测试
最后一个输出显示 'None' ，我该如何删除它？

【解决方案3】：

从数组（列表）的开头删除需要将其后的所有项目向下移动。这意味着以这种方式从列表中从前面开始删除每个元素是一个 O(n^2) 操作。

使用集合可以更有效地做到这一点：

def primes_sieve(limit):
    limitn = limit+1
    not_prime = set()
    primes = []

    for i in range(2, limitn):
        if i in not_prime:
            continue

        for f in range(i*2, limitn, i):
            not_prime.add(f)

        primes.append(i)

    return primes

print primes_sieve(1000000)

... 或者，避免重新排列列表：

def primes_sieve(limit):
    limitn = limit+1
    not_prime = [False] * limitn
    primes = []

    for i in range(2, limitn):
        if not_prime[i]:
            continue
        for f in xrange(i*2, limitn, i):
            not_prime[f] = True

        primes.append(i)

    return primes

【讨论】：

请参阅下面的@Piet Delport 答案以获得优化：将上面的i*2 替换为i*i。

【解决方案4】：

更快：

import time
def get_primes(n):
  m = n+1
  #numbers = [True for i in range(m)]
  numbers = [True] * m #EDIT: faster
  for i in range(2, int(n**0.5 + 1)):
    if numbers[i]:
      for j in range(i*i, m, i):
        numbers[j] = False
  primes = []
  for i in range(2, m):
    if numbers[i]:
      primes.append(i)
  return primes

start = time.time()
primes = get_primes(10000)
print(time.time() - start)
print(get_primes(100))

【讨论】：

【解决方案5】：

我意识到这并不能真正回答如何快速生成素数的问题，但也许有些人会发现这个替代方案很有趣：因为 python 通过生成器提供惰性求值，所以可以完全按照所述实现 Eratosthenes 筛子：

def intsfrom(n):
    while True:
        yield n
        n += 1

def sieve(ilist):
    p = next(ilist)
    yield p
    for q in sieve(n for n in ilist if n%p != 0):
        yield q


try:
    for p in sieve(intsfrom(2)):
        print p,

    print ''
except RuntimeError as e:
    print e

try 块在那里是因为算法一直运行直到它破坏堆栈并且没有尝试阻止显示回溯，将您想要在屏幕外看到的实际输出推送出去。

【讨论】：

不，不是the sieve of Eratosthenes，而是试除法的筛子。即使这样也不是最理想的，因为它不是postponed：任何候选数只需要通过不高于其平方根的素数进行测试。沿着上面链接的答案（后一个）底部的伪代码行实现这一点将为您的代码提供巨大的加速（甚至在您切换到正确的筛子之前）和/因为它将大大减少堆栈的使用 - 所以毕竟你可能不需要你的try 块。
... 另请参阅：more discussion about the "sqrt" issue and its effects、actual Python code for a postponed trial division 和 some related Scala。 --- 如果你自己想出了那个代码，那么对你表示敬意！ :)
有趣，虽然我还不明白为什么我放的东西与 Eratosthenes 的筛子不同。我认为它被描述为将 2 中的所有整数排成一行，然后反复将行中的第一个整数作为素数并删除所有倍数。 “n for n in ilist if n%p != 0”位应该代表剔除倍数。诚然，绝对不是最理想的
n for n in ilist if n%p != 0 测试每个数字 n 在可被p 整除的范围内；但range(p*p, N, p) 直接生成倍数，完全由自身生成，无需测试所有这些数字。

【解决方案6】：

通过结合众多爱好者（包括来自上述 cmets 的 Glenn Maynard 和 MrHIDEn）的贡献，我在 python 2 中提出了以下代码：

def simpleSieve(sieveSize):
    #creating Sieve.
    sieve = [True] * (sieveSize+1)
    # 0 and 1 are not considered prime.
    sieve[0] = False
    sieve[1] = False
    for i in xrange(2,int(math.sqrt(sieveSize))+1):
        if sieve[i] == False:
            continue
        for pointer in xrange(i**2, sieveSize+1, i):
            sieve[pointer] = False
    # Sieve is left with prime numbers == True
    primes = []
    for i in xrange(sieveSize+1):
        if sieve[i] == True:
            primes.append(i)
    return primes

sieveSize = input()
primes = simpleSieve(sieveSize)

在我的机器上以 10 的幂计算不同输入所需的时间是：

3 : 0.3 毫秒
4：2.4 毫秒
5 : 23 毫秒
6 : 0.26 秒
7 : 3.1 秒
8 : 33 秒

【讨论】：

不再需要与 True 或 False 进行比较，因为它们已经是布尔值，
@Copperfield 谢谢！它有助于将速度提高 10-20%。
这个sieve = [True] * (sieveSize+1) 比我的解决方案快，但是sieve[0]/[1] 和xrange(sieveSize+1) at primes[] 并没有改善任何东西。 xrange(2, sieveSize+1) 很好。 :)。也可以代替for i in xrange(2,int(math.sqrt(sieveSize))+1): 使用for i in xrange(2, int((sieveSize+1)**0.5): 好的代码。 :)

【解决方案7】：

使用一点 numpy，我可以在 2 秒多一点的时间内找到所有低于 1 亿的素数。

有两个关键特性需要注意

只为i 删除i 的倍数，直到n 的根
使用x[2*i::i] = False 将i 的倍数设置为False 比显式python for 循环快得多。

这两个显着加快了您的代码速度。对于低于 100 万的限制，没有可察觉的运行时间。

import numpy as np

def primes(n):
    x = np.ones((n+1,), dtype=np.bool)
    x[0] = False
    x[1] = False
    for i in range(2, int(n**0.5)+1):
        if x[i]:
            x[2*i::i] = False

    primes = np.where(x == True)[0]
    return primes

print(len(primes(100_000_000)))

【讨论】：

【解决方案8】：

一个简单的速度技巧：定义变量“primes”时，将步长设置为 2 以自动跳过所有偶数，并将起点设置为 1。

然后你可以进一步优化，而不是 for i in primes，使用 for i in primes[:round(len(primes) ** 0.5)]。这将大大提高性能。此外，您可以消除以 5 结尾的数字以进一步提高速度。

【讨论】：

【解决方案9】：

我的实现：

import math
n = 100
marked = {}
for i in range(2, int(math.sqrt(n))):
    if not marked.get(i):
        for x in range(i * i, n, i):
            marked[x] = True

for i in range(2, n):
    if not marked.get(i):
        print i

【讨论】：

我刚刚测试了你的代码，我看到 dict 解决方案比 list 解决方案慢 2 倍。

【解决方案10】：

这是一个更节省内存的版本（并且：一个适当的筛子，而不是试除法）。基本上，不是保留所有数字的数组，并删除那些不是素数的数字，而是保留一组计数器 - 一个计数器对应于它发现的每个素数 - 并在假定的素数之前跳过它们。这样，它使用的存储空间与素数的数量成正比，而不是最高的素数。

import itertools

def primes():

    class counter:
        def __init__ (this,  n): this.n, this.current,  this.isVirgin = n, n*n,  True
            # isVirgin means it's never been incremented
        def advancePast (this,  n): # return true if the counter advanced
            if this.current > n:
                if this.isVirgin: raise StopIteration # if this is virgin, then so will be all the subsequent counters.  Don't need to iterate further.
                return False
            this.current += this.n # pre: this.current == n; post: this.current > n.
            this.isVirgin = False # when it's gone, it's gone
            return True

    yield 1
    multiples = []
    for n in itertools.count(2):
        isPrime = True
        for p in (m.advancePast(n) for m in multiples):
            if p: isPrime = False
        if isPrime:
            yield n
            multiples.append (counter (n))

您会注意到primes() 是一个生成器，因此您可以将结果保存在列表中，也可以直接使用它们。这是第一个 n 素数：

import itertools

for k in itertools.islice (primes(),  n):
    print (k)

为了完整起见，这里有一个计时器来衡量性能：

import time

def timer ():
    t,  k = time.process_time(),  10
    for p in primes():
        if p>k:
            print (time.process_time()-t,  " to ",  p,  "\n")
            k *= 10
            if k>100000: return

以防万一，我还写了primes() 作为一个简单的迭代器（使用__iter__ 和__next__），它以几乎相同的速度运行。也让我大吃一惊！

【讨论】：

有趣的想法 - 如果将素数计数器存储在最小堆中，它会提高性能（内部循环将是 O(log num_primes) 而不是 O(num_primes)）
为什么？即使它们成堆，我们仍然必须考虑每一个。
如果你将每个素数存储在由它的下一个值作为键的堆中，你只需要查看下一个值为当前值 n 的素数。最大的素数将沉入堆的底部，并且比较小的素数需要评估的次数要少得多。

【解决方案11】：

因为速度，我更喜欢 NumPy。

import numpy as np

# Find all prime numbers using Sieve of Eratosthenes
def get_primes1(n):
    m = int(np.sqrt(n))
    is_prime = np.ones(n, dtype=bool)
    is_prime[:2] = False  # 0 and 1 are not primes

    for i in range(2, m):
        if is_prime[i] == False:
            continue
        is_prime[i*i::i] = False

    return np.nonzero(is_prime)[0]

# Find all prime numbers using brute-force.
def isprime(n):
    ''' Check if integer n is a prime '''
    n = abs(int(n))  # n is a positive integer
    if n < 2:  # 0 and 1 are not primes
        return False
    if n == 2:  # 2 is the only even prime number
        return True
    if not n & 1:  # all other even numbers are not primes
        return False
    # Range starts with 3 and only needs to go up the square root
    # of n for all odd numbers
    for x in range(3, int(n**0.5)+1, 2):
        if n % x == 0:
            return False
    return True

# To apply a function to a numpy array, one have to vectorize the function
def get_primes2(n):
    vectorized_isprime = np.vectorize(isprime)
    a = np.arange(n)
    return a[vectorized_isprime(a)]

检查输出：

n = 100
print(get_primes1(n))
print(get_primes2(n))    
    [ 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]
    [ 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]

比较埃拉托色尼筛的速度和 Jupyter Notebook 上的蛮力。 Eratosthenes 的筛分法比百万元素的蛮力法快 539 倍。

%timeit get_primes1(1000000)
%timeit get_primes2(1000000)
4.79 ms ± 90.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
2.58 s ± 31.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

【讨论】：

您的内循环内容不应该是（考虑到以前的答案和 cmets）if is_prime[i]: is_prime[i*i::2*i]=False 一行吗？

【解决方案12】：

我认为必须可以简单地使用空列表作为循环的终止条件并想出了这个：

limit = 100
ints = list(range(2, limit))   # Will end up empty

while len(ints) > 0:
    prime = ints[0]
    print prime
    ints.remove(prime)
    i = 2
    multiple = prime * i
    while multiple <= limit:
        if multiple in ints:
            ints.remove(multiple)
        i += 1
        multiple = prime * i

【讨论】：

【解决方案13】：

import math
def sieve(n):
    primes = [True]*n
    primes[0] = False
    primes[1] = False
    for i in range(2,int(math.sqrt(n))+1):
            j = i*i
            while j < n:
                    primes[j] = False
                    j = j+i
    return [x for x in range(n) if primes[x] == True]

【讨论】：

【解决方案14】：

我认为这是用 eratosthenes 方法查找素数的最短代码

def prime(r):
    n = range(2,r)
    while len(n)>0:
        yield n[0]
        n = [x for x in n if x not in range(n[0],r,n[0])]


print(list(prime(r)))

【讨论】：

不过，性能绝对是可怕的。它在每次迭代时创建一个全新的列表。

【解决方案15】：

我能想到的最快实现：

isprime = [True]*N
isprime[0] = isprime[1] = False
for i in range(4, N, 2):
    isprime[i] = False
for i in range(3, N, 2):
    if isprime[i]:
        for j in range(i*i, N, 2*i):
            isprime[j] = False

【讨论】：

【解决方案16】：

我刚想出这个。它可能不是最快的，但除了直接添加和比较之外，我没有使用任何其他方法。当然，在这里阻止你的是递归限制。

def nondivsby2():
    j = 1
    while True:
        j += 2
        yield j

def nondivsbyk(k, nondivs):
    j = 0
    for i in nondivs:
        while j < i:
            j += k
        if j > i:
            yield i

def primes():
    nd = nondivsby2()
    while True:
        p = next(nd)
        nd = nondivsbyk(p, nd)
        yield p

def main():
    for p in primes():
        print(p)

【讨论】：

非常好的配方，干净清晰简洁！我会给它加书签。当然，要产生第 100 个素数，nd 链将有 99 层深。但真正需要的只有 10 个。随着我们沿着素数列表走得越远，它变得越来越糟。你能找到解决这个问题的方法吗？ :)
另外，我在这里没有看到任何递归，所以这里也不应该有任何限制。（当然我几乎完全不懂 Python）
一开始我很惊讶，当我得到RecursionError: maximum recursion depth exceeded 异常时。但后来我觉得这有点道理。因为我们可以将生成器视为具有__next__ 函数的对象。所以每个nondivsbyk生成器都是同一个类的对象（只是初始化不同）。假设我们称它为class_nondivsbyk，那么当一个调用另一个时，它实际上是class_nondivsbyk.__next__ 在另一个对象上调用另一个class_nondivsbyk.__next__。
关于第 100 个素数只需要前 10 个素数，所以首先我可以说（只要我不想给出上限）我们需要“收集”方式，因此创建这些生成器似乎是必要的。我现在真的不知道我是否可以“跳过”那些不相关的计算。因为现在，我让每个检查它是否是一个除法器，但如果我把它们放在一边，当数字增加时我需要别的东西来“触发”它们，我不知道如何将它整合到递归中。我还做了一个“平面”版本，我可以在那里看看。谢谢@WillNess
分配nd = nondivsbyk(p, nd)之后的两个nds应该是两个不同的对象。即，nd 是首先引用对象的变量；然后通过函数调用构造新的生成器对象，并将其分配给相同的变量（忘记其旧值）。但在内部，新的生成器对象指的是旧的 - 不同的 - 对象。但正如我所说，我不知道 Python。关于 10 个质数与 100 个质数 - 这里有一个提示：希望每次调用 primes() 都会创建一个单独的新生成器对象。（或者正确的术语是什么？）

【解决方案17】：

我制作了埃拉托色尼筛的单衬版

sieve = lambda j: [print(x) for x in filter(lambda n: 0 not in map(lambda i: n % i, range(2, n)) and (n!=1)&(n!=0), range(j + 1))]

就性能而言，我很确定这无论如何都不是最快的，就可读性/遵循 PEP8 而言，这非常糟糕，但它比任何东西都更新颖。

p>

编辑：请注意，这只是打印筛子并且不会返回（如果您尝试打印它，您将获得一个无列表，如果您想返回，请将列表理解中的 print(x) 更改为“ x"。

【讨论】：

【解决方案18】：

不确定我的代码是否有效，有人愿意评论吗？

from math import isqrt

def isPrime(n):
    if n >= 2: # cheating the 2, is 2 even prime?
        for i in range(3, int(n / 2 + 1),2): # dont waste time with even numbers
            if n % i == 0:
                return False
    return True

def primesTo(n): 
    x = [2] if n >= 2 else [] # cheat the only even prime
    if n >= 2:
        for i in range(3, n + 1,2): # dont waste time with even numbers
            if isPrime(i):
                x.append(i)  
    return x

def primes2(n): # trying to do this using set methods and the "Sieve of Eratosthenes"
    base = {2} # again cheating the 2
    base.update(set(range(3, n + 1, 2))) # build the base of odd numbers
    for i in range(3, isqrt(n) + 1, 2): # apply the sieve
        base.difference_update(set(range(2 * i, n + 1 , i)))
    return list(base)

print(primesTo(10000)) # 2 different methods for comparison
print(primes2(10000))

【讨论】：

【解决方案19】：

获得主要数字的最快方法可能如下：

import sympy
list(sympy.primerange(lower, upper+1))

如果您不需要存储它们，只需使用上面的代码而不转换为list。 sympy.primerange 是一个生成器，所以不会消耗内存。

【讨论】：

请在您的答案正文中解释为什么这是必要的以及它带来了哪些改进使其看起来是一个有意义的答案。

【解决方案20】：

使用递归和海象运算符：

def prime_factors(n):
    for i in range(2, int(n ** 0.5) + 1):
        if (q_r := divmod(n, i))[1] == 0:
            return [i] + factor_list(q_r[0])
    return [n]

【讨论】：