素数筛法算法及其原理

引言

本文介绍部分素数筛法的步骤以及原理，并附带 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