算法-生成质数

前言

获取质数是常见的需求，假如我们想获取 $2 \sim n$（假设 $n \le N$，$N$ 为已知的常数）的所有质数，应该怎么办呢？

两种选择

检查质数

首先不难想到，如果我们想达到这个目的，只需要检查这个范围内的每个数是不是质数就行了。

质数的定义是，除 $1$ 和它本身以外没有其它因数的数，如果我们这么选择，那么在检查一个数 $x$ 时，需要至少判断 $2 \sim \sqrt x$ 中的每个数是否都不能被 $x$ 整除（之所以是 $\sqrt x$ 是因为：如果 $y \gt \sqrt x$ 能被 $x$整除，那么 $\dfrac{x}{y} \lt \sqrt x$ 也能被 $x$ 整除，所以检查小于 $\sqrt x$的部分足矣）。

那么对于 $2 \sim n$ 整体，就要检查 $\displaystyle \sum \sqrt n$ 次。估算时间复杂度约为 $O(n \cdot \sqrt n)$。

筛去合数

有没有更好的做法呢？

已知合数的定义和质数相反，只要存在其它因子就算合数。

那么如果检查每个数是否为合数，就比检查质数方便得多，因为只需要找到一个反例即可。

怎么找合数呢？

埃拉托斯特尼筛法（埃筛）

初步想法

首先对于每个数，将它乘以 $x \in \mathbb N$ 倍得到的一定是合数。

所以我们可以先从 $2$ 开始，得到 $\textcolor{red} 4, \textcolor{blue} 6, 8, \cdots$，再从 $3$ 开始，得到 $\textcolor{blue} 6, 9, 12, \cdots$，由于 $\textcolor{red} 4$ 被筛去，所以跳到 $5 \rightarrow 10, 15, 20, \cdots$，由于 $\textcolor{blue} 6$ 被筛去，所以跳到 $7\cdots$

以C++为例，代码如下（假设长度为 $N$ 的数组不会导致过高的内存占用，$n$ 在int范围内）：

bool is_prime[N + 1];

void gen_primes(int n) {
    for(int i = 2; i <= n; i++) {
        is_prime[i] = true;
    }
    for(int i = 2; i <= n; i++) {
        if(!is_prime[i]) continue;
        for(int j = 2; j <= n / i; j++) {
            is_prime[i * j] = false;
        }
    }
}

两重优化

你肯定注意到这种方式有一定缺陷，比如 $\textcolor{blue} 6$ 同时被 $2 \times 3$ 和 $3 \times 2$ 筛去了，所以可以限制 $j \le i$，这样就不会出现 $3 \times 2$ 这种情况了。

另外，$i \gt \sqrt n$ 时，$i \times j$ 不可能 $ \lt n$，所以只需考虑 $i \in \left[2,\sqrt n \right]$ 即可。

恭喜你发现了埃筛！

bool is_prime[N + 1];
vector<int> primes;

void gen_primes(int n) {
    for(int i = 2; i <= n; i++) {
        is_prime[i] = true;
    }

    const int limit = sqrt(n);
    for(int i = 2; i <= limit; i++) {
        if(!is_prime[i]) continue;
        for(int j = i; j <= n / i; j++) {
            is_prime[i * j] = false;
        }
    }

    for(int i = 2; i <= n; i++) {
        if(is_prime[i]) primes.push_back(i);
    }
}

计算时间复杂度

对于质数 $2$，我们筛去了 $\dfrac n 2$ 个合数，$3$ 筛去了 $\dfrac n 3$ 个合数，$\cdots$

那么总操作次数就是 $\displaystyle \sum_p \dfrac n p $，由于 $\displaystyle \sum \dfrac 1 n \approx \log n$，质数的数量也是 $ \dfrac 1 \log$ 级别，所以总时间复杂度为 $O(n \log \log n)$（注：严谨证明需借助积分，此处仅作直观理解故略去）。

欧拉筛（线性筛）

经过优化之后，我们还是可以发现一些问题。

比如 $30 = 2 \times 15 = 5 \times 6$，就会被筛去两次。

有没有什么办法，能让每个合数都只被筛去一次呢？

筛合数的办法

我们可以将一个合数描述为：最小质因子 $\times$ 其它。

这样就可以对所有的 $i \in \left[2,n\right]$ 乘上质数 $p$，得到所有合数。

提前终止

首先，类似于刚才提到的 $j \le i$，你想到了可以限制 $p \le i$，防止 $15$ 被 $3 \times 5$ 和 $5 \times 3$ 同时筛去。

那么现在产生了一个新问题，$12=6 \times 2 = 4 \times 3$，要怎么避免它同时被 $4$ 和 $6$ 筛去呢？

之所以会有这种情况，是因为 $4$ 有比 $3$ 更小的质因数：$2$，所以才会产生 $4 \times 3 = 2 \times (2 \times 3) = 2 \times 6$ 的情况。

那么既然对于任意一个数 $x$，都有 $4 \times x = 2 \times 2x$，那么 $4x$ 这个数就不需要被 $4$ 或 $x$ 筛去了。

再比如，$18 = 6 \times 3 = 9 \times 2$，根据刚才的推理，因为 $6$ 有更小的质因数 $2$，所以 $18$ 不该被 $6$ 筛去。

所以，如果 $p$ 是 $i$ 的因数，那么就应该停止 比 $p$ 大的质数与 $i$ 相乘（因为 $i$ 有更小的质因数 $p$）。

综上，$p$ 的终止条件有两个（至少一个成立就终止）：

• $p \gt i$
• $i \mod p = 0$

由于当 $i = p$ 时，$i \mod p = 0$，所以只需考虑第二个条件即可。

以下为C++实现

bool is_prime[N + 1];
vector<int> primes;

void gen_primes(int n) {
    for(int i = 2; i <= n; i++) {
        is_prime[i] = true;
    }
    for(int i = 2; i <= n; i++) {
        if(is_prime[i])
            primes.push_back(i);
        for(int p : primes) {
            if(1ll * i * p > n) break;
            is_prime[i * p] = false;
            if(i % p == 0) break;
        }
    }
}

恭喜，你发明了欧拉筛！

这种筛法对每个合数都只筛去了一次，所以时间复杂度是 $O(n)$！

总结

大多数时候，使用埃筛就足够了，如果要追求极致速度，可以选择线性筛。

需要注意，线性筛筛去合数的过程并不是按大小筛的，下面给出一部分筛合数的过程：

合数	$i$	$p$
$4$	$2$	$2$
$6$	$3$	$2$
$9$	$3$	$3$
$8$	$4$	$2$
$10$	$5$	$2$
$15$	$5$	$3$
$25$	$5$	$5$
$12$	$6$	$2$
$14$	$7$	$2$
$21$	$7$	$3$
$16$	$8$	$2$
$18$	$9$	$2$
$27$	$9$	$3$
$20$	$10$	$2$
$22$	$11$	$2$
$24$	$12$	$2$
$26$	$13$	$2$
$28$	$14$	$2$
$30$	$15$	$2$
$\cdots$	$\cdots$	$\cdots$