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 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++代碼顯得更加簡潔緊湊了.

參考