---

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'''

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

优化后的埃氏筛

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
'''

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

这里的优化主要体现在遍历数目中了, 但是由此带来的一个问题就是不能通过一次遍历找出所有的素数, 需要额外遍历.

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
'''

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

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;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

P.S. 这样看反而是C++代码显得更加简洁紧凑了.

参考

---

1. 筛法 - OI Wiki (oi-wiki.org); ↩︎