埃拉托斯特尼筛法 - 寻找素数 Python
- 2024-12-16 08:35:00
- admin 原创
- 149
问题描述：
需要澄清的是，这不是家庭作业问题：）
我想为我正在构建的数学应用程序寻找素数，并遇到了埃拉托斯特尼筛选法。
我已经用 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)
解决方案 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。）
解决方案 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)
解决方案 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
解决方案 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：
利用一点点numpy，我可以在两秒多一点的时间内找到所有小于 1 亿的素数。
有两个关键特征需要注意
i只剪掉 的倍数，i直到 的根ni设置使用False的倍数x[2*i::i] = False比显式的 Python for 循环要快得多。
这两个方法可以显著加快代码速度。对于低于一百万的限制，运行时间几乎感觉不到。
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)))
解决方案 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 块，是因为算法会一直运行直到它耗尽堆栈，而如果没有 try 块，则会显示回溯，从而将您想要看到的实际输出推送到屏幕外。
解决方案 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毫秒
4：2.4毫秒
5：23毫秒
6：0.26秒
7：3.1秒
8 ：33 秒
解决方案 8：
一个简单的速度技巧：当你定义变量“primes”时，将步骤设置为 2 以自动跳过所有偶数，并将起点设置为 1。
然后，您可以进一步优化，使用 for i in primes[:round(len(primes) ** 0.5)] 代替 for i in primes。这将大大提高性能。此外，您可以消除以 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
解决方案 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, "
")
k *= 10
if k>100000: return
如果您想知道的话，我也写了primes()一个简单的迭代器（使用__iter__和__next__），它的运行速度几乎相同。这也让我很惊讶！
解决方案 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)
解决方案 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：
我认为这是用埃拉托斯特尼方法寻找素数的最短代码
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)
解决方案 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”。
解决方案 18：
埃拉托斯特尼筛法各种方法的实证分析与可视化
我在《Python 数据结构与算法》第10章（映射、哈希表和跳过列表）的末尾发现了这个算法。我最终写了三个版本：
这是一个简单的实现，具有令人讨厌的立方运行时间（以及许多其他问题，后来得到解决）。
经过彻底的研究后得到了次优版本，但我仍然想摆脱基于平方的倍数。
性能方面与所选答案相当的快速高效版本，但具有加法倍数（对于那些不太擅长数学的人来说更直观）。
幼稚的 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 为间隔，从范围 100 中的素数到范围 4500 中的素数，进行了 45 次迭代（这是可视化的最佳点，因为尽管图形形状一致，但简单实现的增长使其他三个的可见性相形见绌）。您可以在GitHub gist上调整可视化代码，但这里有一个示例输出：
解决方案 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"
Primes 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"
Enter 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]
解决方案 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))
解决方案 21：
获得主要数字的最快方法可能如下：
import sympy
list(sympy.primerange(lower, upper+1))
如果您不需要存储它们，只需使用上面的代码而无需转换为list。sympy.primerange是一个生成器，因此它不会消耗内存。
解决方案 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]
解决方案 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):,}
'
f'Greatest: {prime_list[-1]:,}
'
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))
解决方案 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))
解决方案 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.6 秒，比我的单核快2.1倍，比您的原始代码快27倍。
换句话说，与你的代码相比，我的多核代码的速度提高了27primes_limit_M = 8倍，真的很多！而且这只适用于 2 核笔记本电脑。如果你有 8 个或更多核心，那么你将获得更大的加速。但请记住，多核机器上的真正加速只会针对相当大的素数限制，尝试 64 百万限制或更大，为此修改我代码中的这一行以设置百万的数量。
只是为了深入讨论细节。我的单核版本几乎和你的代码一样，但使用了Numpy，这使得任何数组操作都非常快，而不是使用带有列表的纯 Python 循环。
多核版本也使用 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()
解决方案 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))
解决方案 27：
如果您希望速度更快，那么如果您有 Nvidia 处理器，您也可以使用 numba 和 cuda。请根据需要随意优化。我们使用 2**24，在 Nvidia 3070 上以 239 毫秒的速度处理约 1600 万个数字。
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)
相关推荐
热门文章
项目管理软件有哪些？
- 2025年20款好用的项目管理软件推荐，项目管理提效的20个工具和技巧
- 2024年开源项目管理软件有哪些？推荐5款好用的项目管理工具
- 2024年常用的项目管理软件有哪些？推荐这10款国内外好用的项目管理工具
- 项目管理软件有哪些？推荐7款超好用的项目管理工具
- 项目管理软件有哪些最好用？推荐6款好用的项目管理工具
- 项目管理软件哪个最好用？盘点推荐5款好用的项目管理工具
- 项目管理软件排行榜：2024年项目经理必备5款开源项目管理软件汇总
- 项目管理必备:盘点2024年13款好用的项目管理软件
- 项目管理软件有哪些，盘点推荐国内外超好用的7款项目管理工具
- 2024项目管理软件排行榜（10类常用的项目管理工具全推荐）
热门标签
云禅道AD
){kind=link}