素数筛法算法及其原理

引言

本文介绍部分素数筛法的步骤以及原理,并附带 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)2 > 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 倍数的过程中,我们是从 p2 开始,而不是从 2·p 开始
  • 之前厄拉多塞筛法的原理中提到过,任意合数都能被不大于 √n 的素数标记。对于 2·p, 3·p, ..., (p-1)·p 这些数,因为这些数开根号之后的结果均小于 p,故这些数如果是合数,就能够被小于 p 的素数标记
  • 故当开始筛 p 的倍数时,小于 p2 的合数都已经被之前的素数筛过了

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 = n1·n2 = (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 = p1·p2···pq,其中 p1 为最小的素数
  • 由于任意小于 p1 的素数都不能整除 p2···pq,所以 n 会在 i = p2···pq 时被筛去
  • 再证明每个合数只会被筛去 1 次:
  • 设合数 n 即被 p1 筛去,也被 p1' 筛去。那么有 n = q·p1 = q'·p1',其中 p1 和 p1' 均是 n 的素因子
  • 不妨设 p1 < p1',故 p1 和 p1' 互素,故有 p1 | q'
  • i = q' 时,遍历素数到 p1 时有 i % p1 == 0,此时跳出循环,不会再遍历到后面的 p1'
  • 故 n 不会被 p1' 筛去,只会被其最小的素因子 p1 筛去

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, k2·seg2]。当 (k+1)·seg ≤ k2·seg2,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 ≤ k2,该条件当且仅当 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 的素数整除。对小于 i2 的数,其要么被小于 i 的素数筛去,要么提前跳出了循环。故将 mp 初始值设为 i2 是合理的,还能够减少运算量

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

参考资料

posted @ 2021-03-27 19:16  kentle  阅读(1912)  评论(0编辑  收藏  举报