素數篩法算法及其原理

引言

本文介紹部分素數篩法的步驟以及原理，並附帶 python 算法的實現
本文介紹的篩法有：

• 厄拉多塞篩法（Eratosthenes Sieve）
• Sundaram 篩法
• 歐拉篩法（Euler Sieve）
• 分段篩法（Segmented Sieve）
• 增量篩（Incremental sieve）
• Atkin 篩法
  
  

厄拉多塞篩法（Sieve of Eratosthenes）

1. 厄拉多塞篩法步驟

給定一個數 n，從 2 開始依次將 \(\sqrt{n}\) 以內的素數的倍數標記為合數
標記完成后，剩余未被標記的數為素數（從 2 開始）
算法原理如下：

1. 讀取輸入的數 n，將 2 至 n 所有整數記錄在表中
2. 從 2 開始，划去表中所有 2 的倍數
3. 由小到大尋找表中下一個未被划去的整數，再划去表中所有該整數的倍數
4. 重復第(3)步，直到找到的整數大於 √n 為止
5. 表中所有未被划去的整數均為素數

2. 厄拉多塞篩法原理

• 首先，先證明這種方法能夠標記所有 2~n 之間的合數。
• 由整數的唯一分解定理知，任意一個大於1的正整數都可以被分解成有限個素數的乘積。因此，任意一個合數都可以看做是一個素數的倍數的形式。
• 對於任意一個合數 n，存在 p, q ，使得 n = p·q （不妨設 p 為素數）
• 同時可有 min(p, q) ≤ p ≤ √n，否則會有 p · q ≥ min(p, q)² ＞ n，矛盾
• 故可知，任意合數都能被不大於 √n 的素數 p 標記
• 其次，顯然，該標記方法並不會錯誤地將素數標記成合數
• 故該方法能且僅能標記所有 2~n 之間的合數，所以剩下未被標記的數就是素數（從 2 開始）

3. 厄拉多塞篩法代碼

from math import sqrt

def sieve_of_eratosthenes(n: int):
    is_prime = [True for _ in range(n + 1)]
    for i in range(2, int(sqrt(n)) + 1):
        if is_prime[i]:
            for j in range(i * i, n + 1, i):
                is_prime[j] = False  # 篩去j
    # 返回從2開始未被篩去的整數
    return [i for i in range(2, n + 1) if is_prime[i]]

• 在篩去素數 p 倍數的過程中，我們是從 p² 開始，而不是從 2·p 開始
• 之前厄拉多塞篩法的原理中提到過，任意合數都能被不大於 √n 的素數標記。對於 2·p, 3·p, ..., (p-1)·p 這些數，因為這些數開根號之后的結果均小於 p，故這些數如果是合數，就能夠被小於 p 的素數標記
• 故當開始篩 p 的倍數時，小於 p² 的合數都已經被之前的素數篩過了
  
  

Sundaram 篩法（Sieve of Sundaram）

1. Sundaram 篩法步驟

Sundaram 篩法使用了一種特殊的方法，繞過所有偶數，篩去所有的奇合數

1. 給定一個正整數 n，令 k = \(\left[\frac{n-1}{2} \right]\)（方括號 [ ] 表示下取整）
2. 將所有不超過 k 的形如 \(i+j+2 \cdot i \cdot j, \;\; i,j \in Z^+\) 的整數標記
3. 記剩下的不超過 k 且未被標記的整數集合為 A
4. 如果 n ≥ 2，則返回 \(\{2\} \; \bigcup \; \{ 2 q +1 | \;q \in A, \; q \not= 0\}\)
   否則返回空集

2. Sundaram 篩法原理

• Sundaram 篩法並沒有考慮偶數的情況，它將奇合數全部篩去，剩余的大於 1 的奇數便都是素數
• Sundaram 篩法的標記數組將數組下標 i 映射到 2i+1，即 0 號位對應整數 1，1 號位對應整數 3
• 對於奇合數 q，有 q = n₁·n₂ = (2i + 1)·(2j + 1) = 2i + 2j + 4ij + 1，對應的下標為 (q - 1) / 2 = i + j + 2ij
• 從上述推導過程也可知，下標 i + j + 2ij 對應的整數 q 必定為合數
• 故可以通過篩去對應下標 i + j + 2ij 的整數來篩去奇合數
• 此外，對於正整數 n，不超過 n 的奇數有 k = \(\left[\frac{n-1}{2} \right]\) 個
• 對於不超過 k 的 i + j + 2ij，在遍歷過程中，若將 i 作為外層循環，j 作為內層循環，那么可以通過如下循環遍歷

for i in range(1, k + 1):
    for j in range(1, k + 1):
        if i + j + 2 * i * j <= k:
            # 篩去 i+j+2ij 對應整數

• 不難發現，i 與 j 具有對稱性，若 j 仍從 1 開始循環，那么會產生較多重復。當 i = m ，j = 1, 2, ..., m-1 時，其實等價於 i = 1, 2, ..., m-1, j = 1，因為二者對應的 i + j + 2ij 的值相等，篩去的都是同一個整數，且后者在 i = m 之前已經被遍歷過。故可以讓 j 從 i 開始遍歷，減少重復遍歷的次數。即：

for i in range(1, k + 1):
    for j in range(i, k + 1):
        if i + j + 2 * i * j <= k:
            # 篩去 i+j+2ij 對應整數

• 對於 i 和 j 的范圍，我們還能夠再優化一下
• 當 j = i 的時候，有 i + j + 2ij = 2i · (i + 1)，這是內層遍歷過程中篩去的下標最小值，同時該值也隨着 i 的循環而遞增。一旦該值超過了整數 k，那么之后的 i + j + 2ij 的值均會超過 k。因此，令 2i · (i + 1) ≤ k，可以解出 i 的最大值為 \(\frac{\sqrt{1+2k}-1}{2}\)。借此可以對外層 i 的循環進行優化
• 同時，可以發現，i + j + 2ij = i + (2i + 1) · j，若將 i 視為定值，則上式是以2i · (i + 1) 為首項，以 (2i + 1) 為公差的等差數列。那么根據這點可以對內層循環進行優化

for i in range(1, int((sqrt(1 + 2 * k) - 1) / 2) + 1):
    for j in range(2 * i *(i + 1), k + 1, 2 * i + 1):
        # 篩去 j 對於整數

3. Sundaram 篩法代碼

from math import sqrt

def sieve_of_sundaram(n: int):
    if n < 2:  # 如果n小於2，則不存在素數，返回空列表
        return []
    else:  # 否則進行線性篩素數計算
        # 即使n=2，算法依然能夠返回正確的值
        k = (n - 1) // 2
        is_prime = [True for _ in range(k + 1)]
        for i in range(1, int((sqrt(1 + 2 * k) - 1) / 2) + 1):
            for j in range(2 * i + 2 * i * i, k + 1, 2 * i + 1):
                is_prime[j] = False  # 篩去對應整數
        return [2] + [2 * i + 1 for i in range(1, k + 1) if is_prime[i]]

歐拉篩法

1. 歐拉篩法步驟

算法原理如下：

1. 讀取輸入的數 n，將 2 至 n 所有整數記錄在表中
2. 從 i=2 開始，如果 i 在表中，則將 i 加入至素數列表中
3. 遍歷當前素數列表，篩去 i 素數倍的數
4. 當遍歷到能整除 i 的素數 p 時，篩去 i·p，停止對素數列表的遍歷
5. 重復 2, 3, 4，直到所有不超過 n 的整數都被遍歷過
6. 素數列表中的元素即為所求的不超過 n 的所有素數

2. 歐拉篩法原理

• 首先證明每個合數都會被篩去：
• 若 n 為合數，設其素因子分解為 n = p₁·p₂···p_q，其中 p₁ 為最小的素數
• 由於任意小於 p₁ 的素數都不能整除 p₂···p_q，所以 n 會在 i = p₂···p_q 時被篩去
• 再證明每個合數只會被篩去 1 次：
• 設合數 n 即被 p₁ 篩去，也被 p₁' 篩去。那么有 n = q·p₁ = q'·p₁'，其中 p₁ 和 p₁' 均是 n 的素因子
• 不妨設 p₁ < p₁'，故 p₁ 和 p₁' 互素，故有 p₁ | q'
• i = q' 時，遍歷素數到 p₁ 時有 i % p₁ == 0，此時跳出循環，不會再遍歷到后面的 p₁'
• 故 n 不會被 p₁' 篩去，只會被其最小的素因子 p₁ 篩去

3. 歐拉篩法代碼

# 輸入正整數n
# 返回不超過n的所有素數
def sieve_of_euler(n):
    primes = []
    is_prime = [True for _ in range(n + 1)]
    for i in range(2, n + 1):
        if is_prime[i]:
            primes.append(i)
        for p in primes:
            k = i * p
            if k > n:
                break
            else:
                is_prime[k] = False
            if i % p == 0:
                break
    return primes

分段篩法（Segmented sieve）

1. 分段篩法步驟

分段篩法分區間段篩取素數，所需要的標記數組大小可以自由指定，數組的大小即為分段的長度。分段 seg 的極端情況是 seg = 1，對應逐個判斷素數

1. 給定正整數 n 和 分段 seg（默認為 √n）
2. 將 [0, n] 區域分段，每段長度為 seg，即分成 [0, seg], (seg, 2·seg], ..., (k·seg, n]
3. 使用厄拉多塞篩法求出 [0, seg] 之間的素數，存在 prime 數組中
4. 對於 prime 數組中的素數 p，標記 p 位於 (seg, 2·seg] 內的倍數
5. 將 (seg, 2·seg] 中未被標記的數加入至 prime 數組中
6. 重復 4, 5 的方法，逐個處理每個分段
7. 當最后一個分段 (k · seg, n] 計算完畢后，返回 prime 數組

2. 分段篩法原理

• 在上文敘述厄拉托塞篩法原理的過程中提到，任意一個合數 n 都能被不超過 √n 的素數篩去
• 對於分段 (k·seg, (k+1)·seg]，在處理該分段時我們已經獲取了 [0, k·seg] 之間的素數，能夠篩去的合數范圍為 [0, k²·seg²]。當 (k+1)·seg ≤ k²·seg²，k ≥ 1 時，篩法必定能夠正確運行
• 解得 $seg \geq \frac{1+k}{k^2} + 1 $，即有 \(seg \geq max \left( \frac{1+k}{k^2} \right) = 2\)
• 故當 seg ≥ 2 時，算法正確
• 而 seg = 1 時，要使篩法正確運行，需要的條件為 k+1 ≤ k²，該條件當且僅當 k=1 時不成立。而 k=1 時，對應的區域為 (1, 2]，即使上式不成立，算法依舊能夠正確地將素數 2 加入至素數數組 prime 中
• 故對於任意正整數 seg，都能夠通過該算法篩出不超過 n 的所有素數
• 使用分段篩法能夠將所需的空間壓縮，降低空間復雜度；但由於在篩素數的過程中產生了大量重復篩數的情況，再加上每個分段篩數過程中對鏈表擴展的操作也很耗時，導致時間復雜度增加。需要根據具體情況考慮分段的長度

3. 分段篩法代碼

from math import sqrt

def segmented_sieve(n: int, seg=None):
    if seg is None:  # seg的默認值
        seg = int(sqrt(n))
    elif seg > n or seg <= 0:  # seg輸入不合法
        seg = n

    # 使用厄拉多塞篩法求出[0,seg]內的素數，具體篩法代碼見上文
    primes = eratosthenes(seg)  
    low, high = seg, seg * 2
    # 循環篩去 (low, high] 對應的合數
    while low < n:
        if high > n:
            high = n
        mark = [True for _ in range(high-low+1)]  # 構建標記列表
        for p in primes:
            start = (low // p + 1) * p  # 從大於low且能被p整除的第一個數開始
            for i in range(start, high+1, p):
                mark[i - low] = False
        # 將未被標記的整數添加至素數列表
        primes.extend([i for i in range(low+1, high+1) if mark[i-low]])
        low, high = low + seg, high + seg
    return primes

Incremental sieve

1. Incremental Sieve 步驟

1. 給定正整數 n
2. 從 2 到 n 遍歷 i，對於每個 i，判斷 i 是否是不超過 √i 的素數的倍數：如果是，則 i 為合數；否則為素數

2. Incremental Sieve 原理

• 不難知道，合數能夠被素數整除。判斷一個正整數 n 是否是合數，可以判斷其是否能被不大於 √n 的數整除。
• 但是除法需要消耗的時間太長了，Incremental Sieve 使用了一種較為巧妙的方法。
• 其額外設置了一個列表 mp，用於存放素數的倍數，mp[k] 表示第 k 個素數（primes[k]）的倍數。每次判斷 n 是否是合數時，比較 n 與每個 mp 的大小：
• 如果 mp[k] < n，則令 mp[k] = mp[k] + primes[k]，直到 mp[k] ≥ n 為止。如果存在 mp[k] = n，則表明 n 是第 k 個素數的倍數，是合數。如果所有的 mp[k] 都不等於 n，則可以說明 n 是素數。
• 采用上述方法可以降低除法的次數，減少判斷所需要的時間。比起厄拉托塞篩法，該方法需要的額外空間大大減少，但缺點是效率較慢。

3. Incremental Sieve 代碼

from math import sqrt

def incremental_sieve(n):
    primes = []  # 存放素數
    mp = []  # 存放素數的倍數
    for i in range(2, n + 1):
        flag = True
        limit = sqrt(i)
        for k in range(len(primes)):
            if primes[k] > limit:
                break
            while mp[k] < i:
                mp[k] = mp[k] + primes[k]
            if mp[k] == i:
                flag = False
                break
        if flag:  # 是素數
            primes.append(i)
            mp.append(i * i)
    return primes

• 在上面代碼中，每次獲取到新素數時，primes 是插入素數，而 mp 則是插入素數的平方。這是由於任一合數 n 都能被不超過 √n 的素數整除。對小於 i² 的數，其要么被小於 i 的素數篩去，要么提前跳出了循環。故將 mp 初始值設為 i² 是合理的，還能夠減少運算量
  
  

Sieve of Atkin

1. Sieve of Atkin 原理

這我真不會。。。先挖個坑，不懂以后會不會記得補上。。。
感興趣的可以去看一看后面關於 Atkin 篩法的參考資料

2. Sieve of Atkin 代碼

def sieve_of_atkin(n):
    if n < 2:
        return []
    elif n == 2:
        return [2]
    elif n == 3:
        return [2, 3]
    else:
        limit = int(sqrt(n))
        is_prime = [False for _ in range(n + 1)]
        for x in range(1, limit + 1):
            for y in range(1, limit + 1):
                # k=4x^2+y^2
                k = 4 * x * x + y * y
                if k <= n and (k % 12 == 1 or k % 12 == 5):
                    is_prime[k] = True ^ is_prime[k]
                # k=3x^2+y^2
                k = 3 * x * x + y * y
                if k <= n and k % 12 == 7:
                    is_prime[k] = True ^ is_prime[k]
                # k=3x^2-y^2
                k = 3 * x * x - y * y
                if x > y and k <= n and k % 12 == 11:
                    is_prime[k] = True ^ is_prime[k]
        for i in range(5, limit + 1):
            if is_prime[i]:
                for j in range(i * i, n + 1, i * i):
                    is_prime[j] = False
        return [2, 3] + [i for i in range(5, n + 1) if is_prime[i]]

參考資料

• 厄拉多塞篩法，歐拉篩法和Segmented Sieve：https://oi-wiki.org/math/sieve/
• 厄拉多塞篩法：https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
• 歐拉篩法：https://web.archive.org/web/20180220153640/http://debug18.com/posts/introduction-to-sieve-method/
• Sundaram 篩法：https://en.wikipedia.org/wiki/Sieve_of_Sundaram
• Segmented Sieve：https://www.geeksforgeeks.org/segmented-sieve/
• Incremental Sieve：https://www.codevamping.com/2019/01/incremental-sieve-of-eratosthenes/
• Atkin 篩法：https://en.wikipedia.org/wiki/Sieve_of_Atkin
• Atkin 篩法：https://fylux.github.io/2017/03/16/Sieve-Of-Atkin/
• Atkin 篩法：https://www.ams.org/journals/mcom/2004-73-246/S0025-5718-03-01501-1/S0025-5718-03-01501-1.pdf