初等数论 —— 筛法
筛法（埃拉托斯特尼筛法、线性筛）
筛法就像是一个过滤器，对筛法算法输入\(n\)，可以过滤出\(1\sim n\)范围内的所有素数，并且大部分的筛法都是作用域\([1,n]\)的，其中的左边界固定为\(1\)或者\(2\)，具体为什么，我们需要接下来了解一下筛法的优化方向才能了解。
暴力筛法
- 针对于小于等于\(n\)有多少个素数的问题，直接的想法就是对于\(1\sim n\)中的每一个数字都进行素性检测。
- 对其中某个数字\(x\)，最直接的素性检测算法就是遍历\([2,\lfloor\sqrt{x}\rfloor]\)，判断\(x\)是否存在除了\(1\)与\(x\)以外的因子，此时的时间复杂度为\(O(\sum_{k=3}^n\sqrt{k})≈O(n^{3\over2})\)。但是具体的，我们还需要看素性检测的算法的时间复杂度。
- 其中存在一些简单的优化方法，例如我们只对奇数进行素性检测……
埃拉托斯特尼筛法
- 为了对暴力筛法进行优化，考虑如果\(x\in[2,{n\over2}]\)，那么\((x,n]\)范围内的\(2x,3x,\cdots\)一定为合数。
- 按从小到大的顺序遍历，然后将当前的数字的倍数都标记为合数，那么运行结束时未被标记的数就是素数。
- 设\(x\in[2,n]\)，若\(x\)是一个素数，\(x\)只有\(1\)和\(x\)两个因子，\(D=[2,x-1]\)遍历的过程中，并不存在\(i\in D\)使得\(i\times j=x\)，素数会被留下；若\(x\)是一个合数，设其可以分解为\(x=small \cdot big\)，则\(small\in [2,\lfloor\sqrt{x}\rfloor]\subseteq D\)，\(x\)一定会被\(small\)标记过了，即非素数会被筛去。
- 我们在遍历到\(x\)的时候，如果\(x\)不是素数，那么在其\(small\)被遍历的过程中，我们已经对其倍数（\(small×big,2×small×big,\cdots\)）进行了标记，因此直接跳过；如果\(x\)是素数，是不是需要逐个标记\(2x,3x\cdots\)为非素数？答案是并不需要，我们只需要标记\(x^2,(x+1)\cdot x\cdots\)，因为\(2x\)在遍历到\(2\)的时候给筛掉了，\(3x,\cdots,(x-1)\cdot x\)同理。
vector<int> primes;
vector<bool> isPrime(n + 1, true);
void Eratosthenes(int n) {
isPrime[0] = isPrime[1] = false;
for (int i = 2; i <= n; i++) {
if (isPrime[i]) {
primes.push_back(i);
// 防止 i * i 溢出 int 最大正数为负数
if ((long long)i * i > n) {
continue;
}
for (int j = i * i; j <= n; j += i) {
isPrime[j] = false;
}
}
}
}
- 以上就是\(Eratosthenes\)筛法，时间复杂度为\(O(n\ log\ log\ n)\)
- 上面我们可以看到，进行\(x^2\)的筛选前，我们会判定其是否超出了\(n\)的范围，但因为\(x^2\)是单调递增的，因此当第一次出现这种状况之后，我们就不需要对之后的数据进行判定是否进行筛选了，即只需要对\([2,\lceil\sqrt{n}\rceil]\)执行筛法：
vector<int> primes;
vector<bool> isPrime(n + 1, true);
void Eratosthenes(int n) {
isPrime[0] = isPrime[1] = false;
for (int i = 2; i * i <= n; i++) {
if (isPrime[i]) {
for (int j = i * i; j <= n; j += i) {
isPrime[j] = false;
}
}
}
for (int i = 2; i <= n; i++) [
if (isPrime[i]) {
primes.push_back(i);
}
]
}
- 我们可以只对奇数进行筛选做进一步优化。
- 位操作减少内存占用思想
分块筛法
- 按照上面的思路，我们发现在“筛”的过程中，我们只需要执行到\(\lceil\sqrt{n}\rceil\)，不需要一直保留整个
isPrime[1..n]，为了进行筛选，只需要保留到\(\sqrt{n}\)的素数就足够了，即primes[1..sqrt(n)]，并将整个范围分成块，每个块分别进行筛选。 - 分块筛法的核心思想就是缓存友好性。当\(n\)很大（如\(10^9\)）时，开辟一个大小为\(n\)的布尔数组会导致严重的内存超出。设\(s\)是一个常数，它决定了块的大小，那么我们就有了\(\lceil{n\over s}\rceil\)个块，而块\(k(k=0...\lfloor{n\over s}\rfloor)\)包含了区间\([ks,ks+s-1]\)中的数字。我们依次处理块，也就是说，对于每个块\(k\)，我们将遍历所有素数（从\(1\)到\(\sqrt{n}\)）并使用它们进行筛选。
int countPrime(int n) {
constexpr static int S = 10000;
// 预处理基础素数 [1, sqrt(n)]
vector<int> primes;
int nsqrt = sqrt(n);
vector<bool> isPrime(nsqrt + 1, true);
for (int i = 2; i <= nsqrt; i++) {
if (isPrime[i]) {
primes.push_back(i);
for (int j = i * i; j <= nsqrt; j += i) isPrime[j] = false;
}
}
int result = 0;
vector<bool> block(S);
// 初始化分块边界 [start, start + S - 1]
// block[j] = start + j
for (int k = 0; k * S <= n; k++) {
fill(block.begin(), block.end(), true); // 每次处理新块，假设它们都是素数
int start = k * S;
// 将全局素数映射到局部块
for (int p : primes) {
// 寻找当前块中，素数 p 的最小倍数，求一最小 m ，使得 m × p ≥ start
// m = ⌈start / p⌉
int startIndex = (start + p - 1) / p;
int j = max(startIndex, p) * p - start; // 最少从 p * p 开始筛， - start 是映射到 block 内
for (; j < S; j += p) block[j] = false;
}
if (k == 0) block[0] = block[1] = false;
for (int i = 0; i < S && start + i <= n; i++) {
if (block[i]) result++;
}
}
return result;
}
#include <iostream>
#include <vector>
#include <cmath>
int countPrime(unsigned int L, unsigned int R) {
int S = R - L + 1;
int Rsqrt = std::sqrt(R);
std::vector<int> primes;
std::vector<bool> isPrime(Rsqrt + 1, true);
for (int i = 2; i <= Rsqrt; i++) {
if (isPrime[i]) primes.push_back(i);
for (int j = i * i; j <= Rsqrt; j += i) isPrime[j] = false;
}
int result = 0;
std::vector<bool> block(S, true);
for (int p : primes) {
int startIndex = (L + p - 1) / p;
int j = std::max(startIndex, p) * p - L;
for (; j < S; j += p) block[j] = false;
}
if (L == 1) block[0] = false;
for (int i = 0; i < S; i++) {
if (block[i]) result++;
}
return result;
}
int main() {
unsigned int L = 0;
unsigned int R = 0;
scanf("%d %d", &L, &R);
int ans = countPrime(L, R);
printf("%d\n", ans);
return 0;
}
线性筛
- 在埃拉托斯特尼筛法中，我们通过避免对\(2x,3x,\cdots\)的重复判定来优化算法。但它还不是最优的，其还会存在多余的操作，以上面的代码为例：
n = 30
i = 2: Mark 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30
i = 3: Mark 9, 12, 15, 18, 21, 24, 27, 30
i = 5: Mark 25, 30
- 其中\(12\)被\(2\)和\(3\)各标记\(1\)次，\(30\)被\(2\)、\(3\)、\(5\)各标记\(1\)次。线性筛的优化思路就是让每个合数仅被它的最小的质因子标记\(1\)次，即\(12\)仅被\(2\)标记。
- 由算数基本定理，每一个合数\(C\)都可以唯一地分解为若干素数的成绩。如果我们强制每一个合数都只能被其最小的素因子筛掉，那么每个数字就只会被访问一次，从而达到严格的\(O(n)\)。
- 设当前的遍历数为\(i\)，已知集合\(P=\{p_0,p_1,p_2,\cdots\}\)，我们要做的事情是，用\(i\)去乘以比它小的素数\(p_j\)，构造出一个合数\(C=i\times p_j\)，从而将其筛掉。由于我们是从小到大遍历素数\(p_j\)的，在触发
break之前，\(p_j\)一定是\(C\)的最小素因子。 - 假设\(i\pmod{p_j}\equiv 0\)，这意味着\(i\)可以分解为\(i=k\times p_j\)，如果我们不
break，继续往下执行，循环会取到下一个更大的素数\(p_{j+1}\)，此时我们会标记合数\(C_{next}=i\times p_{j+1}\)，我们将\(i\)展开，代入\(C_{next}\)的表达式：\[C_{next}=(k\times p_j)\times p_{j+1}=(k\times p_{j+1})\times p_j \]显然，\(C_{next}\)的最小素因子是\(p_j\)，而不是\(p_{j+1}\)，如果这里用\(p_{j+1}\)把\(C_{next}\)筛掉，这就违反了“只用最少素因子筛选”的原则。这个\(C_{next}\)应该在未来当外部循环到\(i'=k\times p_{j+1}\)时，再乘以素数\(p_j\)被筛掉。
vector<int> primes;
vector<bool> isPrime(n + 1, true);
void Euler(int n) {
for (int i = 2; i <= n; i++) {
if (isPrime[i]) {
primes.push_back(i);
}
for (int j = 0; j < primes.size() && i * primes[j] <= n; j++) {
isPrime[i * primes[j]] = false;
if (i % primes[j] == 0) {
// i % primes[j] == 0
// 意味着 i 在 primes[j] 就已经被筛掉了
// 这还意味着 i 的倍数 k * i = k * primes[j] * xxx 也被筛掉了
break;
}
}
}
}
- 线性筛的时间复杂度为\(O(n)\).
欧拉筛求素数的同时也得到了每个数的最小质因子，即当
break发生时，\(i\)的最小素因子就找到了。
- 参考资料：OI-Wiki
浙公网安备 33010602011771号