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

Question

只是为了澄清，这不是一个家庭作业问题:)

我想为我正在构建的数学应用程序找到素数遇到了埃拉托斯特尼筛法方法。

我已经用 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)

Answer 1

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

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

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

正确的算法（使用列表而不是字典）类似于：（

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) 处开始非素数标记的算法优化，而不是它的两倍。）

Answer 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)

Answer 3

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

您可以使用 set: ... 更有效地完成此操作

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

Answer 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))

Answer 5

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

有两个关键特征需要注意，

• 仅针对 i 删除 i 的倍数，直到 n 的根号
• 设置 i< 的倍数使用 x[2*i::i] = False 将 /code> 转换为 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)))

Answer 6

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

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 块在那里是因为该算法运行直到它耗尽堆栈并且没有
尝试阻止显示回溯，将您想要在屏幕上看到的实际输出推开。

Answer 7

通过结合许多爱好者的贡献（包括上述评论中的 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 ms
• 4 : 2.4 毫秒
• 5 : 23 毫秒
• 6 : 0.26 秒
• 7 : 3.1 秒
• 8 : 33 秒

Answer 8

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

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

Answer 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

Answer 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__），并且它以几乎相同的速度运行。也让我很惊讶！

Answer 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 上比较埃拉托斯特尼筛法和暴力破解的速度。埃拉托色尼筛法比暴力法快 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)

Answer 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

Answer 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]

Answer 14

我认为这是用埃拉托色尼方法查找素数的最短代码

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

Answer 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

Answer 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)

Answer 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 而言，这非常糟糕，但它更新颖的长度比什么都重要。

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

Answer 18

埃拉托斯特尼筛法的各种方法的实证分析和可视化

我在第 10^th 章（映射、哈希表和跳过列表）末尾发现了该算法
“数据结构和算法” Python”。我最终编写了三个版本：

• 一个幼稚的实现，带有令人讨厌的立方运行时（以及许多其他问题，后来解决了）。
• 经过一些 彻底的研究，但我仍然想摆脱基于平方的倍数。
• 一个快速高效的版本，在性能方面与选定的答案相当，但具有加法倍数（对于那些不太了解的人来说更直观）数学倾向）。

Naive SOE

在所有问题中，返回语句中的额外空间使用是最糟糕的。

def naive_sieve(m: int):
    BA = [True] * m
    for i, k in zip(range(2, m + 1), range(len(BA))):
        if BA[k] is False: continue
        for j in range(2, i):
            if i % j == 0:
                BA[k] = False
                f = k + j
                while f < len(BA):
                    BA[f] = False
                    f += j
                break
    return [i for i,j in zip(range(2, m + 1), BA) if j is True]

次优 SOE

有了巨大的改进，但仍然缺乏空间使用，并且内循环间隔进展不友好。

def suboptimal_sieve(m: int):
    BA = [True] * m
    for i, k in zip(range(2, m + 1), range(2, len(BA))):
        if BA[k] is False: continue
        for j in range(i**2, m, i):
            BA[j] = False
    return [i for i,j in zip(range(2, m + 1), BA[2:]) if j is True]

快速 SOE

易于理解（注意小数指数的使用与根的数学定义明确一致）
和 @Pi Delport 的答案一样高效。

def fast_sieve(m: int):
    BA = [True] * m
    rtm = int(m**(1/2)) + 1
    for i in range(2, len(BA)):
        if BA[i]:
            yield i
            if i < rtm:
                f = i
                while f < len(BA):
                    BA[f] = False
                    f += i

实证分析

要比较所有三种实现以及 @Pi Delport 所选的答案，
我从 100 范围内的素数到 4500 范围内的素数，以 100 为间隔，进行了 45 次迭代
（这是可视化的最佳点，因为尽管图形形状保持一致，
天真的实现的增长使其他三个的可见性相形见绌）。
您可以在 GitHub gist 上调整可视化代码，但这里有一个示例输出：

Answer 19

Bitarray 和 6k±1 的大小和速度

5 秒内 10 亿 3 毫秒内 2^24

我使用 bitarray 来减少占用空间。 bitarray 还允许这样的语句：

 p[4::2] = False  # Clear multiples of 2

筛子代码将在我的旧 i7 上大约 5 秒内完成 10 亿大小的筛子

$ python prime_count_test.py 
Creating new sieve of 100,000,000 primes.
 Make sieve for 100,000,000 took 0:00:00.306957 
Count primes to: >1000000000
Creating new sieve of 1,000,000,000 primes.
 Make sieve for 1,000,000,000 took 0:00:04.734377 
From 1 to 1,000,000,000 there are 50,847,534  primes
 Search took 0:00:04.856540 

Primes from 1,000,000 to 1,000,200 are 
1000003, 1000033, 1000037, 1000039, 1000081, 1000099, 1000117, 1000121, 1000133, 1000151, 1000159, 1000171, 1000183, 1000187, 1000193, 1000199

Enter a number < 1,000,000,000  to test for prime> 10017
10,017  is not prime

这是我的筛子：

from datetime import timedelta
import time
from bitarray import bitarray

def get_primes(a_prime, start=0, end=100):
    """Return list of prime numbers in the range [start, end] using a_prime bitarray."""
    return [i for i in range(start, end + 1) if a_prime[i]]

def make_sieve(size):
    """Create a sieve of Eratosthenes up to the given size."""
    sieve_time = time.time()
    print(f"Creating new sieve of {size:,} primes.")
    limit = int(1 + size**0.5) + 2
    p = bitarray(size)  # One bit per value
    p.setall(True)
    p[0] = p[1] = False
    p[4::2] = False  # Clear multiples of 2
    p[9::3] = False  # Clear multiples of 3

    for i in range(5, limit, 6):  # Process only numbers of the form 6k-1 and 6k+1
        h = i + 2  # 6k+1
        if p[i]:  # If 6k-1 is prime
            p[3 * i::2 * i] = False  # Clear multiples of 6k-1
        if p[h]:  # If 6k+1 is prime
            p[3 * h::2 * h] = False  # Clear multiples of 6k+1
            
    print(f"Make sieve for {len(p):,} took {str(timedelta(seconds=time.time() - sieve_time))}")
    return p

sieve_size = 10**8
p = make_sieve(sieve_size)

# Input for counting primes up to n
n = int(input("Count primes to: >"))
start_time = time.time()

if n <= sieve_size:
    prime_count = p[1:n + 1].count()
else:
    p = make_sieve(n)
    prime_count = p[1:n + 1].count()
    sieve_size = len(p)

print(f"From 1 to {n:,} there are {prime_count:,} primes")
print(f"Search took {str(timedelta(seconds=time.time() - start_time))}")

# Get 200 primes from the sieve starting at 1,000,000
print(f"\nPrimes from 1,000,000 to 1,000,200 are:")
print(*get_primes(p, 10**6, 10**6 + 200), sep=', ')

# Test if a specific number is prime
test_p = int(input(f"\nEnter a number < {sieve_size:,} to test for prime> "))
print(f"{test_p:,} {'is prime' if p[test_p] else 'is not prime'}")

仅返回素数：

from datetime import timedelta
import time
from bitarray import bitarray

def get_primes(a_prime, start=0, end=100):
    """Return list of prime numbers in the range [start, end] using a_prime bitarray."""
    return [i for i in range(start, end + 1) if a_prime[i]]

def make_sieve(size):
    """Create a sieve of Eratosthenes up to the given size."""
    sieve_time = time.time()
    print(f"Creating new sieve of {size:,} primes.")
    limit = int(1 + size**0.5) + 2
    p = bitarray(size)  # One bit per value
    p.setall(True)
    p[0] = p[1] = False
    p[4::2] = False  # Clear multiples of 2
    p[9::3] = False  # Clear multiples of 3

    for i in range(5, limit, 6):  # Process only numbers of the form 6k-1 and 6k+1
        h = i + 2  # 6k+1
        if p[i]:  # If 6k-1 is prime
            p[i * i::2 * i] = False  # Clear multiples of 6k-1
        if p[h]:  # If 6k+1 is prime
            p[h * h::2 * h] = False  # Clear multiples of 6k+1
            
    print(f"Make sieve for {len(p):,} took {str(timedelta(seconds=time.time() - sieve_time))}")
    
    return [i for i in range(len(p)) if p[i]]
    
start_time = time.time()
n = 2**24
p = make_sieve(n)
print(f"From 1 to {n:,} there are {len(p):,} primes")

print(p[0:200])

输出 2**24 搜索大约需要 3 毫秒：

$ python prime_count_test.py 
Creating new sieve of 16,777,216 primes.
Make sieve for 16,777,216 took 0:00:00.034416
From 1 to 16,777,216 there are 1,077,871 primes
[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, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223]

Answer 20

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

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

Answer 21

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

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

如果您不需要存储它们，只需使用上面的代码，无需转换为列表。 sympy.primerange 是一个生成器，因此它不消耗内存。

Answer 22

使用递归和海象运算符：

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]

Answer 23

基本筛选的

使用 numpy 进行 速度非常快。可能是最快的实现

# record: sieve 1_000_000_000 in 6.9s (core i7 - 2.6Ghz)
def sieve_22max_naive(bound):
    sieve = np.ones(bound, dtype=bool)  # default all prime
    sieve[:2] = False  # 0, 1 is not prime

    sqrt_bound = math.ceil(math.sqrt(bound))

    for i in range(2, sqrt_bound):
        if sieve[i]:
            inc = i if i == 2 else 2 * i
            sieve[i * i:bound:inc] = False

    return np.arange(bound)[sieve]


if __name__ == '__main__':
    start = time.time()
    prime_list = sieve_22max_naive(1_000_000_000)
    print(f'Count: {len(prime_list):,}\n'
          f'Greatest: {prime_list[-1]:,}\n'
          f'Elapsed: %.3f' % (time.time() - start))

分段筛（使用更少的内存）

# find prime in range [from..N), base on primes in range [2..from)
def sieve_era_part(primes, nfrom, n):
    sieve_part = np.ones(n - nfrom, dtype=bool)  # default all prime

    limit = math.ceil(math.sqrt(n))
    # [2,3,5,7,11...p] can find primes < (p+2)^2
    if primes[-1] < limit - 2:
        print(f'Not enough base primes to find up to {n:,}')
        return

    for p in primes:
        if p >= limit: break

        mul = p * p
        inc = p * (2 if p > 2 else 1)
        if mul < nfrom:
            mul = math.ceil(nfrom / p) * p
            (mul := mul + p) if p > 2 and (mul & 1) == 0 else ...  # odd, not even

        sieve_part[mul - nfrom::inc] = False

    return np.arange(nfrom, n)[sieve_part]
    # return np.where(sieve_part)[0] + nfrom
    # return [i + nfrom for i, is_p in enumerate(sieve_part) if is_p]
    # return [i for i in range(max(nfrom, 2), n) if sieve_part[i - nfrom]]


# find nth prime number, use less memory,
# extend bound to SEG_SIZE each loop
# record: 50_847_534 nth prime in 6.78s, core i7 - 9850H 2.6GHhz
def nth_prime(n):
    # find prime up to bound
    bound = 500_000
    primes = sieve_22max_naive(bound)

    SEG_SIZE = int(50e6)

    while len(primes) < n:
        # sieve for next segment
        new_primes = sieve_era_part(primes, bound, bound + SEG_SIZE)
        # extend primes
        bound += SEG_SIZE
        primes = np.append(primes, new_primes)

    return primes[n - 1]


if __name__ == '__main__':
    start = time.time()
    prime = nth_prime(50_847_534)
    print(f'{prime:,} Time %.6f' % (time.time() - start))

Answer 24

这是我的解决方案，与维基百科相同

import math


def sieve_of_eratosthenes(n):
    a = [i for i in range(2, n+1)]
    clone_a = a[:]
    b = [i for i in range(2, int(math.sqrt(n))+1)]
    for i in b:
        if i in a:
            c = [pow(i, 2)+(j*i) for j in range(0, n+1)]
            for j in c:
                if j in clone_a:
                    clone_a.remove(j)
    return clone_a


if __name__ == '__main__':
    print(sieve_of_eratosthenes(23))

Answer 25

感谢您提出有趣的问题！

现在我从头开始写了两个版本的经典埃拉托斯特尼筛法。

一种是单核（CPU 核心），另一种是多核（使用所有 CPU 核心）。

两个（单核和多核）版本的主要加速是由于使用了官方 Python 包 Numpy。通过命令python -m pip install numpy安装一次。多核版本的另一个巨大的加速是由于使用了所有CPU核心，与单核版本相比，N个核心的加速几乎是N倍。

现在我只有两核笔记本电脑。我的程序产生了以下计时：

Computing primes less than 8M
Number of CPU cores 2
Original time 72.57 sec
Simple time 5.694 sec
Simple boost (compared to original) 12.745x
Multi core time 2.642 sec
Multi core boost (compared to simple) 2.155x

上面的 Original 表示您问题正文中的代码，这是您已经优化的第二个代码。 Simple 是我的单核版本。 Multi 是我的多核版本。

正如您在上面看到的，在原始代码上计算小于 800 万的素数花费了 72 秒，我的单核版本花费了 5.7 秒，这比您的代码快12.7 倍倍，而我的 2 核版本花费了 5.7 秒。 2.6 秒，比我的单核快 2.1 倍 倍，比您的原始代码快 27 倍 倍。

换句话说，与您的代码相比，我的多核代码的速度提高了 27 倍 倍，真的很多！这仅适用于 2 核笔记本电脑。如果您有 8 个或更多核心，那么您将获得更大的加速。但请记住，多核机器上的真正加速仅适用于相当大的素数限制，请尝试 6400 万限制或更大，为此修改代码中的 primes_limit_M = 8 行以设置百万数量。

只是为了详细说明。我的单核版本几乎与您的代码类似，但使用 Numpy ，这使得任何数组操作都非常快，而不是使用纯带列表的 pythonic 循环。

多核版本也使用 Numpy，但将筛选范围的数组分割成与计算机上的 CPU 核心数量相同的块，每块数组具有相同的大小。然后，每个 CPU 核心在其自己的数组部分中设置布尔标志。此技术仅在达到内存 (RAM) 速度限制之前提供加速，因此在 CPU 核心数量增长的某个点之后，您不会获得额外的加速。

默认情况下，我在多核版本中使用所有 CPU 核心，但您可以通过设置比计算机上拥有的更少的核心来进行实验，这可能会带来更好的加速，因为不能保证大多数核心都会给出最快的结果。通过将 cpu_cores = mp.cpu_count() 行更改为 cpu_cores = 3 等内容来调整内核数量。

在线尝试！

def SieveOfEratosthenes_Original(end):
    import numpy as np
    limit = end - 1
    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 np.array([i for i in primes if primes[i]==True], dtype = np.uint32)

def SieveOfEratosthenes_Simple(end):
    # https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
    import numpy as np
    composites = np.zeros((end,), dtype = np.uint8)
    composites[:2] = 1
    for p, c in enumerate(composites):
        if c:
            continue
        composites[p * p :: p] = 1
    return np.array([p for p, c in enumerate(composites) if not c],
        dtype = np.uint32)

def SieveOfEratosthenes_Task(small_primes, begin, end):
    import numpy as np
    composites = np.zeros((end - begin,), dtype = np.uint8)
    offsets = np.full((len(small_primes),), begin, dtype = np.uint32)
    offsets = small_primes - offsets % small_primes
    offsets[offsets == small_primes] = 0
    for off, p in zip(offsets, small_primes):
        composites[off :: p] = 1
    return np.array([begin + i for i, c in enumerate(composites) if not c],
        dtype = np.uint32)

def SieveOfEratosthenes_MultiCore(end, *, nthreads = None):
    import math, multiprocessing as mp, numpy as np
    end_small = math.ceil(math.sqrt(end)) + 1
    small_primes = SieveOfEratosthenes_Simple(end_small)
    if nthreads is None:
        nthreads = mp.cpu_count()
    block = (end - end_small + nthreads - 1) // nthreads
    with mp.Pool(nthreads) as pool:
        return np.concatenate([small_primes] + pool.starmap(SieveOfEratosthenes_Task, [
            (small_primes, min(end_small + ithr * block, end), min(end_small + (ithr + 1) * block, end))
            for ithr in range(nthreads)]))

def Test():
    import time, numpy as np, multiprocessing as mp
    primes_limit_M = 8
    cpu_cores = mp.cpu_count()
    end = primes_limit_M * 2 ** 20
    print(f'Computing primes less than {primes_limit_M}M')
    print('Number of CPU cores', cpu_cores, flush = True)
    
    tim_orig = time.time()
    res_orig = SieveOfEratosthenes_Original(end)
    tim_orig = time.time() - tim_orig
    print('Original time', round(tim_orig, 3), 'sec', flush = True)
    
    tim_simple = time.time()
    res_simple = SieveOfEratosthenes_Simple(end)
    tim_simple = time.time() - tim_simple
    print('Simple time', round(tim_simple, 3), 'sec', flush = True)
    
    assert np.all(res_orig == res_simple)
    print(f'Simple boost (compared to original) {tim_orig / tim_simple:.3f}x')
    
    tim_multi = time.time()
    res_multi = SieveOfEratosthenes_MultiCore(end, nthreads = cpu_cores)
    tim_multi = time.time() - tim_multi
    print('Multi core time', round(tim_multi, 3), 'sec', flush = True)
    
    assert np.all(res_simple == res_multi)
    print(f'Multi core boost (compared to simple) {tim_simple / tim_multi:.3f}x')
    
if __name__ == '__main__':
    Test()

Answer 26

import math

def atkin_sieve(limit):
   primes = [False] * (limit + 1)
   square_limit = int(math.sqrt(limit))

# Отметить все числа, делящиеся нацело на 2, 3 или 5
for i in range(1, square_limit + 1):
    for j in range(1, square_limit + 1):
        num = 4 * i**2 + j**2
        if num <= limit and (num % 12 == 1 or num % 12 == 5):
            primes[num] = not primes[num]

        num = 3 * i**2 + j**2
        if num <= limit and num % 12 == 7:
            primes[num] = not primes[num]

        num = 3 * i**2 - j**2
        if i > j and num <= limit and num % 12 == 11:
            primes[num] = not primes[num]

# Удалить кратные квадратам простых чисел
for i in range(5, square_limit):
    if primes[i]:
        for j in range(i**2, limit + 1, i**2):
            primes[j] = False

# Вернуть список простых чисел
return [2, 3] + [i for i in range(5, limit) if primes[i]]

# Пример использования
print(atkin_sieve(100))

Answer 27

如果您寻求更快的速度，您也可以使用 numba 和 cuda
你有一个 Nvidia 处理器。请随意根据需要进行优化。我们使用 2**24，在 Nvidia 3070 上，大约 1600 万个数字，耗时 239 毫秒。

from numba import cuda
import numpy as np

@cuda.jit
def findprimes(primes, sqrt_limit):
    index = cuda.grid(1)
    if index >= primes.size or index < 2:
        return
    if index <= sqrt_limit:
        if primes[index]:
            for multiple in range(index*index, primes.size, index):
                primes[multiple] = False

def fast_sieve_gpu(m):
    primes = np.ones(m, dtype=bool)
    primes[0] = primes[1] = False
    sqrt_limit = int(np.sqrt(m))
    d_primes = cuda.to_device(primes)
    threads_per_block = 128
    blocks_per_grid = (m + threads_per_block - 1) // threads_per_block
    findprimes[blocks_per_grid, threads_per_block](d_primes, sqrt_limit)
    primes = d_primes.copy_to_host()
    prime_numbers = np.nonzero(primes)[0]
    return prime_numbers

m = 2**24 # 16777216
prime_numbers = fast_sieve_gpu(m)
print(prime_numbers)
%timeit fast_sieve_gpu(m)

# Output (this is 2**24 which is 16x of 2**20) : 
# [       2        3        5 ... 16777183 16777199 16777213]
# 239 ms ± 2.66 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)