202606221117素数筛选方法

素数筛选方法

什么是素数？

大于 1 的自然数，如果只能被 1 和该数本身整除，不能被其他自然数整除，则称为素数（质数）。相对的，大于 1 的自然数若不是素数，则称之为合数。

因此，若 $a$ 是素数，则

• $a\in N^{+},a>1$；
• 1 既不是素数，也不是合数；
• 素数是在自然数即 $N^{+}$ 的范畴内讨论的。

性质

算数基本定理

每个大于 1 的整数均可以携程一个以上的素数之乘积，且除了质因数的排序不同外是唯一的。

例如，

$$\begin{array}{rl}23244& =2\cdot 2\cdot 3\cdot 13\cdot 149\\ & =2^2\cdot 3\cdot 13\cdot 149\end{array}$$

可见，素数可被认为是构建自然数大厦的“砖块”。

欧几里得定理

存在无穷多个素数。

筛选方法

试除法

如何判断一个数是不是素数，最简单的方法就是根据定义，枚举 2 到 n-1 之间的数，看是否能整除 n，进而判断 n 是否为素数。这种方法叫试除法。

def IsPrime(n):
    if n < 2:
        return False
    for i in range(2, n):
        if n % i == 0:
            return False
    return True

更进一步，对自然数 $n$ 进行因数分解，

$$n=a\cdot b⟹a\in \left[2,\lfloor \sqrt{n}\rfloor \right],b\in \left[\lceil \sqrt{n}\rceil ,n-1\right]$$

也就是说，若 $n$ 可被分解为两个因数之积，一定有较小的因数在 $2$ 到 $\lfloor \sqrt{n}\rfloor$ 之间，较大的因数在 $\lceil \sqrt{n}\rceil$ 到 $n-1$ 之间。因此，我们在尝试找到 $a$ 时，如果 $n$ 真的可被因数分解，那么对应的 $b$ 一定在较大的区间尝试过了。或者说，当 $a$ 跑到第二个区间时，此时的 $b$ 会跑到第一个区间，而第一个区间我们已经尝试过了。所以，我们只需要枚举 $2$ 到 $\sqrt{n}$ 这个区间就行。

def IsPrime(n):
    if n < 2:
        return False
    for i in range(2, int(math.sqrt(n))):
        if n % i == 0:
            return False
    return True

埃拉托斯特尼筛法（sieve of Eratosthenes）

简称埃氏筛。用来筛选出给定自然数 $n$ 以内的素数。其基本步骤是从最小的素数 2 开始，将该素数的所有倍数标记成合数，而下一个尚未被标记的最小自然数 3 即是下一个素数。如此重复这一过程，将各个素数的倍数标记为合数并找出下一个素数，最终便可找出一定范围内所有素数。

N = 5000
is_prime = [True] * N
is_prime[0] = is_prime[1] = False
 
def Eratosthenes(n):
    for i in range(2, n + 1):
        if is_prime[i]:
            # 从 2i 开始，以步长 i 增加，即 2i，3i，4i，...
            for j in range(2 * i, n + 1, i):
                # 是 i 的倍数均不是素数
                is_prime[j] = False
 
# 生成素数表
Eratosthenes(N)
 
def IsPrime(n): # n < N
    return is_prime[n]

上面的筛选过程，有 2 个小优化，

1. j 可以从 i*i 开始，因为 i 的 2 倍到 i - 1 倍，已经筛过了。换句话说，对于 $\{a\cdot i|a=2,3,\dots ,i-1\}$，它也可以写作 $\{i\cdot a|a<i\}$，即小于 i 的数 a 的 i 倍。显然，这个数在处理 a 的时候已经筛选过了。
2. 套用试除法中的结论，当 $i>\sqrt{n}$ 时，$i^2>n$，此时已超出范围，无需再筛。

优化后，

def Eratosthenes(n):
    for i in range(2, n + 1):
        if is_prime[i]:
            if i * i > n:
                break
            for j in range(i * i, n + 1, i):
                is_prime[j] = False

上述代码仍有优化空间，如只筛奇数。

线性筛法

埃氏筛法有重复计算，例如，12 会被 2 标记：4，6，8，10，12，也会被 3 标记：9，12。30 会被

• 2 标记：4, 6, …, 28, 30
• 3 标记：9, 12, 15, …, 27, 30
• 5 标记：25, 30

所以要想办法把这些重复计算给处理掉。这就引出了线性筛法，又称欧拉筛法。

primes = []
is_prime = [True] * N
is_prime[0] = is_prime[1] = False
 
def EulerSeive(n):
    for i in range(2, n + 1):
        if is_prime[i]:
            # 记录目前为止已经筛得的素数
            primes.append(i)
        for p in primes:
            if i * p > n:
                break
            is_prime[i * p] = False
            if i % p == 0:
                """
                i % p == 0
                表明 i 之前被 p 筛过了，p 是 i 的素因子。可以写作 i = k * p.
                而 primes 中的素数是从小到大的，所以 i 乘上后面更大的质数的结果
                
                    i * p2 = k*p * p2 = k*p2 * p = k' * p
 
                一定会被 p 的更大倍数筛掉，
                就不需要在这里先筛一次，所以这里直接 break 掉
                """
                break

设 $primes=[p_1,p_2,p_3,...]$ 是从小到大的素数序列。考虑当 $i\mathrm{mod}{p}_k=0$ 时，对于 $p_x,x>k$，有

$$i\cdot p_{k+1}=(a\cdot p_k)\cdot p_{k+1}=\hat{a}\cdot p_k,\text{where }\hat{a}=a\cdot p_{k+1}$$

这就说明 i 和 $p_k$ 之后素数的乘积，一定会被 $p_k$ 的倍数筛掉，因此这里可以先不筛。

注意，当 i % p == 0 时，p 是 i 的最小素因子。

运行步骤

i=2, p=2, seive 4
i=3, p=2, seive 6
i=3, p=3, seive 9
i=4, p=2, seive 8
i=5, p=2, seive 10
i=5, p=3, seive 15
i=5, p=5, seive 25
i=6, p=2, seive 12
i=7, p=2, seive 14
i=7, p=3, seive 21
i=8, p=2, seive 16
i=9, p=2, seive 18
i=9, p=3, seive 27
i=10, p=2, seive 20
i=11, p=2, seive 22
i=12, p=2, seive 24
i=13, p=2, seive 26
i=14, p=2, seive 28
i=15, p=2, seive 30
i=16, p=2, seive 32
i=17, p=2, seive 34
i=18, p=2, seive 36
i=19, p=2, seive 38
i=20, p=2, seive 40
[2, 3, 5, 7, 11, 13, 17, 19]

References

1. 维基百科 - 素数
2. 埃拉托斯特尼筛法
3. 素数筛法 - OI Wiki