#Python 中这段非常简短的代码尝试模拟前 N 个自然数的“埃拉托斯特尼筛法”，并具有 (0) 脚本简短的约束； (1) 最小化“if 语句”和“for/while 循环”； (2) CPU 时间方面的效率。
import numpy as np
N = 10**5
a = np.array(range(3,N,2))
for j in range(0, int(round(np.sqrt(N),0))):
a[(a!=a[j]) & (a%a[j] == 0)] = 0
a = a[a!=0]
a = [2]+list(a)
在 Intel Core I5 上，它返回第一个中的质数：
有人愿意在上述限制内分享 CPU 时间更高效的代码吗？
埃拉托色尼的实际 NumPy 筛看起来像这样：
def sieve(n):
flags = numpy.ones(n, dtype=bool)
flags[0] = flags[1] = False
for i in range(2, n):
# We could use a lower upper bound for this loop, but I don't want to bother with
# getting the rounding right on the sqrt handling.
if flags[i]:
flags[i*i::i] = False
return numpy.flatnonzero(flags)
它维护一个“可能是素数”标志的数组，并直接取消设置与素数倍数对应的标志，而不需要测试整除性，特别是对于不能被当前正在处理的素数整除的数字。
您所做的是试除法，您只需遍历并测试数字是否可以被候选除数整除。即使很好地实现了试除，也需要比筛子进行更多的操作，并且操作成本更高。您的实现所做的工作甚至比这更多，因为它考虑非素数候选除数，并且因为它不断对它应该已经知道是素数的数字执行整除性测试。
我决定尝试一下，并创建了一个由 @user2357112 发布的进一步优化的 NumPy 版本，支持 Monica，它使用 Numba JIT 来加快速度。
import numba
import numpy
import timeit
import datetime
@numba.jit(nopython = True, parallel = True, fastmath = True, forceobj = False)
def sieve (n: int) -> numpy.ndarray:
primes = numpy.full(n, True)
primes[0], primes[1] = False, False
for i in numba.prange(2, int(numpy.sqrt(n) + 1)):
if primes[i]:
primes[i*i::i] = False
return numpy.flatnonzero(primes)
if __name__ == "__main__":
timestart = timeit.default_timer()
print(sieve(1000000000))
timestop = timeit.default_timer()
timedelta = (timestop - timestart)
print(f"Time Elapsed: {datetime.timedelta(seconds = timedelta)}")
else:
pass
在我的笔记本电脑上，我在
1e90:00:10.3786860:01:27.0599631.94 秒，10.000.000
def sieve_eratosthene(limit):
primes = [True] * (limit+1)
iter = 0
while iter < limit**0.5 :
if iter < 2:
primes[iter]= False
elif primes[iter]:
for i in range(iter*2, limit+1, iter):
primes[i] = False
iter+=1
return(x for x in range(limit+1) if primes[x])
这是使用 NumPy 执行此操作的简单方法。它有点类似于OP的索引而不是在主循环内再次循环的想法，但没有除法检查，只有切片。
它也类似于user2357112支持Monica答案，但是这个只考虑奇数，这使得它更快。
这是典型的仅奇数筛子：https://stackoverflow.com/a/20248491/8094047.
最后我们有一个 bool 数组，我们可以用它来检查 n 以下的数字是否是质数（偶数除外，您可以在没有数组的情况下检查它们，使用 & 1 左右）
平均。 n = 20000000 的时间：0.1063s
import numpy as np
n = 20000000
isprime = np.ones(n, dtype=np.bool)
# odd only sieve
i = 3
while (i * i < n):
if isprime[i]:
isprime[i * i:n:2 * i] = False
i += 2
# test
print(isprime[97]) # should be True
这使用 Numpy 并跳过切片：
def make_sieve(size):
sieve_time = time.time()
print(f"Creating new sieve of {size:,} primes.")
limit = int(1 + size**0.5)+2
p = np.ones(size,dtype=np.int8)
p[0] = p[1] = 0
s = p[4::2]
s[:] = 0 # clear 4,6,8,10
s = p[9::3]
s[:] = 0 # clear 9,12,15...
for i in range(5,limit,6): # all 6k-1
h = i+2 # 6k+1
if p[i]:
s = p[3*i::2*i] # Get a view
s[:] = 0 # Clear it
if p[h]:
s = p[3*h::2*h] # Get a view
s[:] = 0 # Clear it
print(f" Make sieve for {len(p):,} took {str(timedelta(seconds=time.time()-sieve_time))} ")
return p
我计算素数并查找或显示素数，例如：
p = make_sieve(n)
c = np.sum(p[1:n+1])
print(f"From 1 to {n:,} there are {c:,} primes")
def get_primes(x,start=0,size=200):
return [i for i in range(start,start+size+1) if x[i]]
def isPrime(x,n):
return x[n] == 1
输出如下：
Creating new sieve of 100,000,000 primes.
Make sieve for 100,000,000 took 0:00:00.475777
Count primes to: >1000000000
Creating new sieve of 1,000,000,000 primes.
Make sieve for 1,000,000,000 took 0:00:05.382128
From 1 to 1,000,000,000 there are 50,847,534 primes