写在前面
一直想总结一下素数的筛法, 总是抽不开空, 下面用C++和Python实现, 简单讲一下思路, 主要参考了oi-wiki1, 一个打竞赛的大佬们创建的知识集合.
Eratosthenes筛法
思路很简单, 就是通过遍历, 找出已经是素数的数的所有倍数, 将其标记为合数, 那么一趟全部遍历下来, 就能得到所有的素数了.
from time import time
from numba import jit
n = int(1e6)
@jit(nopython=True)
def Eratosthenes(n):
p = 0 # the number of prime
prime = [] # save prime
is_prime = [True] * (n + 1)
is_prime[0] = is_prime[1] = False
for i in range(2, n + 1):# 从第一个素数2开始
if is_prime[i]:# 如果是素数
prime.append(i)#加入素数列表
p = p + 1# 素数个数+1
if i * i <= n:# 不使数组越界, 这两行代码可以不写, 直接进入while判断
j = i * i
while j <= n:
is_prime[j] = False
j = j + i#这里进行的就是倍数的标记, 通过j+=i方式添加倍数
print(prime, "\nthe number of prime is: ", p)
s = time()
Eratosthenes(n)
e = time()
print(f"Eratosthenes: time is {e-s}s")
'''the number of prime is: 78498
Eratosthenes: time is 0.23691391944885254s'''
优化后的埃氏筛
from time import time
# from numba import jit
n = int(1e6)
# @jit(nopython=True)
def Eratosthenes(n):
ans = []#存放素数
cnt = 0
is_prime = [True] * (n + 1)#标记合数
is_prime[0] = is_prime[1] = False# 初始条件
for i in range(2, int(n**.5) + 1):#优化的部分
"""
这里由于只判断了前sqrt(n)个数(这已经能够标记出所有的合数了),
就只能通过第二次遍历得到的bool数组`is_prime`来找出所有的素数,
而不能如前一种方法通过一次遍历来完成素数的存储/计数
"""
if is_prime[i]:
j = i * i
while j <= n:
is_prime[j] = False
j += i
for j in range(n + 1):
if is_prime[j]:
ans.append(j)
cnt += 1
print(ans, "\nthe number of prime is: ", cnt)
s = time()
Eratosthenes(n)
e = time()
print(f"Eratosthenes: time is {e-s}s")
'''the number of prime is: 78498
Eratosthenes: time is 0.23756694793701172s
with jit:
the number of prime is: 78498
Eratosthenes: time is 0.3765590190887451s
'''
这里的优化主要体现在遍历数目中了, 但是由此带来的一个问题就是不能通过一次遍历找出所有的素数, 需要额外遍历.
Euler筛法(线性筛法)
时间复杂度O(N), 相当于上面的代码的进一步优化, 主要思路还是筛法, 但是设置了终止条件, 使得不会进行重复遍历, 提高了运行效率, 这在C++中有所体现.
但是在使用Python进行测试的时候, 埃氏筛竟然是最快的, 不管用没用numba优化, 都是一样慢, 感觉可能是Python对数组有一定优化, 使用最后面给出的C++代码就不会有问题.
from time import time
# from numba import jit
# @jit(nopython=True)
def euler():
MAXN = int(1e6)
pri = []#存储素数
vis = [True] * (MAXN + 1)#标记合数:False
cnt = 0#计数
for i in range(2, MAXN):
if vis[i]:#如果是素数, 存储并计数
pri.append(i)
cnt += 1
for j in range(cnt):#这个循环需要注意
if i * pri[j] > MAXN:# 判断数组越界
break
vis[i * pri[j]] = False # 倍数标记为合数
if i % pri[j] == 0: # 防止重复标记
# 这步是Euler筛法的核心
"""
可以举这样一个例子: 12=3x4=2x6,在素数列表为[2,3],i=4时已进行标记
所以在i=6时候,i%pri[j]=6%2==0,这时候就不会重复标记12了,
同理可证其他像12这样有多素因子的合数不会被重复标记,这就完成了对埃氏筛的优化
"""
break
print(pri, "\nthe number of prime is: ", cnt)
s = time()
euler()
e = time()
print(f"euler: time is {e-s}s")
'''the number of prime is: 78498
euler: time is 0.5525140762329102s
with numba jit:
the number of prime is: 78498
euler: time is 0.40300703048706055s
'''
C++代码
最后整合一下C++版本的代码, 如下:
#include <iostream>
using namespace std;
constexpr int maxn = 1e8 + 10;
bitset<maxn> pri;
int primes[maxn];
void era() {
int N = 1e8, cnt = 0;
double s = clock();
for (int i = 2; i * i <= N; ++i) {
if (!pri[i]) {
for (int j = i * i; j <= N; j += i)
pri[j] = 1;
}
}
for (int i = 2; i <= N; i++)
if (!pri[i])
cnt++;
double e = clock();
printf("%d\ntime = %.0lftic", cnt, e - s);
/*5761455
time = 4252883tic[Finished in 4.9s]*/
}
void euler() {
int N = 1e8, cnt = 0;
double s = clock();
for (int i = 2; i <= N; ++i) {
if (!pri[i])
primes[++cnt] = i;
for (int j = 1; i * primes[j] <= N; j++) {
pri[i * primes[j]] = 1;
if (i % primes[j] == 0)
break;
}
}
double e = clock();
printf("%d\ntime = %.0lftic", cnt, e - s);
/*5761455
time = 2730051tic[Finished in 3.5s]*/
}
int main(int argc, char const* argv[]) {
// era();
euler();
return 0;
}
P.S. 这样看反而是C++代码显得更加简洁紧凑了.