tags: Algorithm Number-Theory

写在前面

一直想总结一下素数的筛法, 总是抽不开空, 下面用C++和Python实现, 简单讲一下思路, 主要参考了oi-wiki¹, 一个打竞赛的大佬们创建的知识集合.

Eratosthenes筛法

思路很简单, 就是通过遍历, 找出已经是素数的数的所有倍数, 将其标记为合数, 那么一趟全部遍历下来, 就能得到所有的素数了.

from time import time
from numba import jit
n = int(1e6)@jit(nopython=True)
def Eratosthenes(n):p = 0  # the number of primeprime = []  # save primeis_prime = [True] * (n + 1)is_prime[0] = is_prime[1] = Falsefor i in range(2, n + 1):# 从第一个素数2开始if is_prime[i]:# 如果是素数prime.append(i)#加入素数列表p = p + 1# 素数个数+1if i * i <= n:# 不使数组越界, 这两行代码可以不写, 直接进入while判断j = i * iwhile j <= n:is_prime[j] = Falsej = 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 = 0is_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 * iwhile j <= n:is_prime[j] = Falsej += ifor j in range(n + 1):if is_prime[j]:ans.append(j)cnt += 1print(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.23756694793701172swith 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)#标记合数:Falsecnt = 0#计数for i in range(2, MAXN):if vis[i]:#如果是素数, 存储并计数pri.append(i)cnt += 1for j in range(cnt):#这个循环需要注意if i * pri[j] > MAXN:# 判断数组越界breakvis[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这样有多素因子的合数不会被重复标记,这就完成了对埃氏筛的优化"""breakprint(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.5525140762329102swith numba jit:the number of prime is:  78498
euler: time is 0.40300703048706055s
'''

C++代码

最后整合一下C++版本的代码, 如下:

#include 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++代码显得更加简洁紧凑了.

参考

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！