NumPy 06-Python实验六-4 埃拉托斯特尼筛法-Numpy.py
埃拉托斯特尼筛法（Sieve of Eratosthenes，简称埃氏筛）
✅ 核心思想
要找出不大于 N 的全部素数：先假设所有数都是素数，然后不断把已知素数的倍数标记为 “不是素数”，筛完剩下未被标记的就是素数。
素数定义：大于 1，只能被 1 和自身整除的整数。
📜 手工演示（举例：N=10，找出 ≤10 的素数）
数字列表：2, 3, 4, 5, 6, 7, 8, 9, 10
- 取第一个素数 2：划掉所有 2 的倍数 → 4、6、8、10 剩下：
2, 3, 5, 7, 9 - 取下一个没划掉的数 3：划掉 3 的倍数 →9 剩下：
2,3,5,7 - 下一个数是 5。\(5>\sqrt{10}\approx3.16\)，停止循环 最终素数：2,3,5,7
关键结论：外层循环只需要到 \(\boldsymbol{\sqrt{N}}\)。 任何合数x，一定存在一个因子 \(\le\sqrt{x}\)；如果遍历到\(\sqrt{N}\)还没筛掉 x，x 一定是素数。
# 小于200万的素数之和： import time import numpy as np start=time.time() n=2000000 prime=np.arange(n+1,dtype=np.uint64) r=int(np.sqrt(n)) for j in range(2,r+1): if prime[j] != 0: prime[j*j:n+1:j]=0 print(np.int64(prime.sum()-1)) end=time.time() print("time:",end-start) # 142913828922 # time: 0.10690665245056152
代码逐行详细讲解
算法：埃拉托斯特尼筛法（埃氏筛，Sieve of Eratosthenes） 目标：求出200 万以内所有素数的总和，用 NumPy 数组加速筛法，比纯 Python 列表循环快很多。 答案结果：142913828922
import time
import numpy as np
import time：计时，统计程序运行耗时import numpy as np：使用 numpy 数组做筛子，利用向量化赋值提速
start=time.time()
记录程序开始的时间戳（秒）。
n=2000000
上限：我们要找 小于等于 2000000 的所有素数。
prime=np.arange(n+1,dtype=np.uint64)
np.arange(n+1)：生成数组[0,1,2,3,4,...,2000000]dtype=np.uint64：无符号 64 位整数，防止大数求和溢出
这个数组
prime就是筛子： 初始：prime[i] = i后面如果判定i不是素数，就把prime[i]置为0； 最后数组里不为 0 的元素，就是素数。
r=int(np.sqrt(n))
\(\sqrt{2000000}\approx1414.21\)，取整 r=1414。
埃氏筛原理：只需要筛到\(\sqrt{n}\)即可。 如果一个数n存在大于\(\sqrt{n}\)的因子，那必然对应一个小于\(\sqrt{n}\)的因子，所以外层循环只需要到\(\sqrt{n}\)。
for j in range(2,r+1):
if prime[j] != 0:
prime[j*j:n+1:j]=0
外层循环，j从 2 遍历到 1414：
if prime[j] != 0：说明j还没有被置 0 → j 是素数prime[j*j : n+1 : j] = 0- 切片：从 \(j^2\) 开始，到 n 为止，步长 j
- 含义：把 j 的所有倍数（\(j^2,j^2+j,j^2+2j\cdots\)）全部标记为 0（非素数）
为什么从
j*j开始，不从2*j？ 小于\(j^2\)的 j 倍数，早就被更小的素数筛掉了，直接从\(j^2\)起步减少重复工作，优化速度。
✨关键点：
prime[...] = 0是NumPy 向量化批量赋值，不是 Python for 循环逐个赋值，这是代码速度很快的核心原因！
print(np.int64(prime.sum()-1))
prime.sum()：数组全部元素求和。- 数组里
prime[0]=0，prime[1]=1。1 不是素数，但初始数组 prime [1]=1，会被算进总和，所以减去 1。
举例：数组元素
[0,1,2,3,0,5,0,7...]，sum 包含了 1，素数是 2,3,5,7，总和 = sum-1。np.int64：转换类型，防止 uint64 输出格式问题。
end=time.time()
print("time:",end-start)
记录结束时间，计算耗时并打印。
算法整体流程总结
- 创建数组
0~2000000，数组下标代表数字本身； - j 从 2 到√2000000 循环；
- 如果 j 是素数，就把 j² 之后所有 j 倍数全部标记为 0；
- 筛完之后，数组中非 0 的数字就是素数；
- 全部相加，减去多余的 1，得到所有素数之和。
重点细节 & 坑点
- 为什么用
uint64？ 200 万素数总和非常大（1429 亿），普通 32 位 int 存不下，会整数溢出，所以用 64 位整型。 - 筛法起点
j*j优化： 朴素写法会写成从2*j开始，结果正确但是慢；j*j是经典优化。 - NumPy 切片赋值是向量化操作
prime[j*j:n+1:j]=0，底层 C 批量赋值，没有 Python 内部循环，这就是为什么 0.1 秒就能跑完 200 万。纯 Python 列表实现埃氏筛，耗时会高很多。 - 必须
sum()-1prime [1]=1，1 不是素数，但没有被筛掉，求和时需要减掉。prime [0]=0 不影响。
简单手写小例子（n=10，方便理解筛法）
初始prime = [0,1,2,3,4,5,6,7,8,9,10]
r = sqrt(10)=3
j=2，prime[2]≠0：prime[4::2]=0 → 数组变为[0,1,2,3,0,5,0,7,0,9,0]
j=3，prime[3]≠0：prime[9::3]=0 → 数组变为[0,1,2,3,0,5,0,7,0,0,0]
sum = 0+1+2+3+0+5+0+7+0+0+0 =18
sum-1 =17
素数：2,3,5,7 总和=17 ✔
对比：原生 Python 列表版本（速度慢很多）
import time
n=2000000
prime = [1]*(n+1)
prime[0]=prime[1]=0
r=int(n**0.5)
start=time.time()
for j in range(2,r+1):
if prime[j]:
for k in range(j*j,n+1,j):
prime[k]=0
s=sum(i for i, is_p in enumerate(prime) if is_p)
print(s)
print(time.time()-start)
原生版本内层是 Python 循环，耗时通常在 1 秒级别，比 numpy 版本慢十倍左右。
浙公网安备 33010602011771号