筛法（埃筛、线性筛）学习笔记

大家好，欢迎来到IT知识分享网。

一、问题引入

求 $1$ ~ $n$ 闭区间内的所有质数。

二、问题分析

1. 朴素的素数筛法

很自然的一种想法就是对 $1$ 到 $n$ 之间的所有数进行一次素数检验，不过显然这种方式很低效，大致复杂度为 $O(n\sqrt{n})$。

2. Eratosthenes筛法

2.1 Eratosthenes筛法的代码实现

#include<bits/stdc++.h> using namespace std; const int maxn = 1e8; int prime[maxn], cnt; bool vist[maxn + 1]; void Eratothenes(int n){ 
    for (int i = 2; i <= n; i ++){ 
    if (!vist[i]){ 
    prime[++ cnt] = i; for (long long j = 1ll * i * i; j <= n; j += i){ 
    vist[j] = true; } } } printf("%d\n", cnt); } int main(){ 
    Eratothenes(maxn); return 0; } 

以上是Eratothenes筛法的代码，其复杂度为 $O(nloglogn)$。实测运行时间大致如下。

2.2 Eratosthenes筛法的时间效率优化

2.2.1 减小筛法上界

需要求出 $1$ 到 $n$ 的所有质数，并不需要取质数进行筛选，只需要取小于 $\sqrt{n}$ 的质数筛选即可。变化如下：

void Eratothenes(){ 
    int i; for (i = 2; i * i <= n; i ++){ 
    // 循环终止条件变化 if (!vist[i]){ 
    prime[++ cnt] = i; for (long long j = 1ll * i * i; j <= n; j += i){ 
    vist[j] = true; } } } // 统计后续的质数 for (; i <= n; i ++) if (!vist[i]) prime[++ cnt] = i; printf("%d\n", cnt); } 

不过优化并不会减少埃筛的时间复杂度，只是减少运算次数而已。

2.2.2 只筛奇数

由于非 $2$ 的偶数都是合数，所以不用考虑奇数，只关心奇数即可。
这样，相当于少了一半的操作。变化如下：

void Eratothenes(int n){ 
    int i; prime[++ cnt] = 2; for (i = 3; i * i <= n; i += 2){ 
    // 步长变为2 if (!vist[i]){ 
    prime[++ cnt] = i; for (long long j = 1ll * i * i; j <= n; j += 2 * i){ 
    vist[j] = true; } } } for (; i <= n; i += 2) // 步长变为2 if (!vist[i]) prime[++ cnt] = i; printf("%d\n", cnt); } 

加上两种优化后会提高埃筛的时间效率，实测运行时间大致如下。

2.3 Eratosthenes筛法的空间效率优化

2.3.1 质数个数的渐进

由质数个数函数的渐近 $\pi (x)\approx \frac{n}{\ln n}$ 可知，质数数组的长度可以比 $n$ 小，大致是 $\frac{n}{\ln n}$ 量级，使用不定长数组(vector)可以稍稍节约空间。

2.3.2 压位

vist数组的每一位其实只需要一个bit，因此可以采用vector<bool> 或者bitset优化空间，从 $n$ byte优化至 $n$ bit。

2.3.3 分块存储

由于只需要取小于 $\sqrt{n}$ 的质数进行筛选，因此，只需要保存 $1$ 到 $\sqrt{n}$ 的vist数组，然后对后面的数分成 $\lceil \frac{n}{s}\rceil$ 个大小为 $s$ 块，用前几个质数筛每一个块里的数，因此每次只需要 $\lceil \frac{n}{s}\rceil$ 大小的数组，总空间复杂度为 $O(\sqrt{n}+s)$。优化后代码如下：

void Eratothenes(int n){ 
    int i; vector <int> prime; int sqrt_n = sqrt(n); vector <bool> vist(sqrt_n + 1, false); // 筛出1~sqrt(n)的所有质数 for (int i = 2; i <= sqrt_n; i ++){ 
    if (!vist[i]){ 
    prime.push_back(i); for (long long j = i * i; j <= sqrt_n; j += i) vist[j] = true; } } // 分块 const int sz = 1e4; int cnt = 0; for (int num = 0; num * sz <= n; num ++){ 
    vector <bool> blk(sz, false); int l = num * sz; // 区间左端点 if (num == 0) blk[0] = blk[1] = true; // 特判 for (int i = 0; i < prime.size(); i ++){ 
    // 使用每一个质数进行筛选 int p = prime[i]; int st = (l + p - 1) / p; // 块内首个可以被p整除的数 for (int j = max(st * p, p * p) - l; j < sz; j += p) blk[j] = true; } for (int i = 0; i < sz && i + l <= n; i ++) if (!blk[i]) cnt ++; } printf("%d\n", cnt); } 

3. Euler筛法

3.1 Euler筛法的代码实现

#include <bits/stdc++.h> using namespace std; const int maxn = 1e8; int prime[maxn], cnt; bool vist[maxn + 1]; void Euler(int n){ 
    for (int i = 2; i <= n; i ++){ 
    if (!vist[i]) prime[++ cnt] = i; for (int j = 1; j <= cnt && i * prime[j] <= n; j ++){ 
    vist[i * prime[j]] = true; if (i % prime[j] == 0) break; } } printf("%d\n", cnt); } int main(){ 
    Euler(maxn); return 0; } 

以上是Euler筛法的代码，其复杂度为 $O(n)$，实测运行时间略慢于埃筛。这点令我十分惊讶，究其原因可能是因为Euler筛法需要多次取模，导致常数稍大。

3.2 Euler筛法的应用

由于每次通过两数的乘积筛去合数，因此Euler筛法可以用来求某些积性函数的值。

3.2.1 Euler筛法求欧拉函数的值

• case 1: $gcd(i,prime[j])=1$，由积性函数的性质有：
  $$\phi (i\ast prime[j])=\phi (i)\ast \phi (prime[j])=\phi (i)\ast (prime[j]-1)$$
  
• case 2: $gcd(i,prime[j])=prime[j]$，根据 $i\ast prime[j]$ 和 $i$ 的唯一分解、欧拉函数的性质有：
  $$\phi (i\ast prime[j])=(i\ast prime[j])\ast {\Pi }_{i=1}^k(1-\frac{1}{p_i})=prime[j]\ast (i\ast {\Pi }_{i=1}^k(1-\frac{1}{p_i}))=prime[j]\ast \phi (i)$$
  

代码如下。

#include <bits/stdc++.h> using namespace std; const int maxn = 10000; int phi[maxn]; int prime[maxn], cnt; bool vist[maxn + 1]; void Euler(int n){ 
    phi[1] = 1; for (int i = 2; i <= n; i ++){ 
    if (!vist[i]) prime[++ cnt] = i, phi[i] = i - 1; for (int j = 1; j <= cnt && i * prime[j] <= n; j ++){ 
    vist[i * prime[j]] = true; if (i % prime[j] != 0) phi[i * prime[j]] = phi[i] * (prime[j] - 1); else { 
    phi[i * prime[j]] = phi[i] * prime[j]; break; } } } } int main(){ 
    Euler(maxn); for (int i = 2; i <= maxn; i ++) printf("%d: %d\n", i, phi[i]); return 0; } 

3.2.2 Euler筛法求约数个数函数的值

约数个数函数，用 $d(x)$ 表示，其值为 $x$ 的约数个数。根据 $x$ 的唯一分解和乘法原理可得：
若 $x={\Pi }_{a=1}^k{p}_a^{c_a}$，则 $d(x)={\Pi }_{i=a}^k(c_a+1)$。
易证约数个数函数为积性函数，尝试通过Euler筛法计算，分类讨论。

• case 1: $gcd(i,prime[j])=1$，由积性函数性质有：
  $$d(i\ast prime[j])=d(i)\ast d(prime[j])=2d(i)$$
  
• case 2: $gcd(i,prime[j])=prime[j]$，由 $i\ast prime[j]$ 和 $i$ 的唯一分解、约数个数函数的性质有（由于 $prime[j]$ 为 $i$ 的最小质因子，不妨设 $i=(prime[j]{)}^{c_1}\ast {\Pi }_{a=2}^k{p}_a^{c_a}$）：
  $$d(i\ast prime[j])=(c_1+2)\ast {\Pi }_{a=2}^k(c_a+1)=\frac{c_1+2}{c_1+1}\ast {\Pi }_{a=1}^k{c}_a=\frac{c_1+2}{c_1+1}d(i)$$
  

因此，要计算 $d(x)$ 的值，只需要记录最小质因子对应的次数 $c_1$ 后用同样的方式计算即可。
代码如下。

#include <bits/stdc++.h> using namespace std; const int maxn = 1e6; int d[maxn], c[maxn]; int prime[maxn], cnt; bool vist[maxn + 1]; void Euler(int n){ 
    d[1] = 1; for (int i = 2; i <= n; i ++){ 
    if (!vist[i]) prime[++ cnt] = i, d[i] = 2, c[i] = 1; for (int j = 1; j <= cnt && i * prime[j] <= n; j ++){ 
    vist[i * prime[j]] = true; if (i % prime[j] != 0){ 
    c[i * prime[j]] = 1; d[i * prime[j]] = d[i] * 2; } else { 
    c[i * prime[j]] = c[i] + 1; d[i * prime[j]] = d[i] / c[i * prime[j]] * (c[i * prime[j]] + 1); break; } } } } int main(){ 
    Euler(maxn); for (int i = 1; i <= 100; i ++) printf("%d: %d\n", i, d[i]); return 0; } 

3.2.3 Euler筛法求约数和函数的值

约束和函数，记作 $sumd(x)$，其值为 $x$ 的所有约数的和。根据 $x$ 的唯一分解和加乘原理可得：
若 $x={\Pi }_{a=1}^k{p}_a^{c_a}$，则 $sumd(x)={\Pi }_{a=1}^k({\Sigma }_{b=0}^{c_a}{p}_a^b)$。
同样易证，约束和函数为积性函数，尝试通过Euler筛法计算，分类讨论。

• case 1: $gcd(i,prime[j])=1$，由积性函数性质有：
  $$sumd(i\ast prime[j])=sumd(i)\ast sumd(prime[j])=sum(d)\ast (1+prime[j])$$
  
• case 2: $gcd(i,prime[j])=prime[j]$，由 $i\ast prime[j]$ 和 $i$ 的唯一分解、约数个数函数的性质有（由于 $prime[j]$ 为 $i$ 的最小质因子，不妨设 $i=(prime[j]{)}^{c_1}\ast {\Pi }_{a=2}^k{p}_a^{c_a}$）：
  $$sumd(i\ast prime[j])={\Sigma }_{b=0}^{c_1+1}prime[j{]}^b\ast {\Pi }_{a=2}^k({\Sigma }_{b=0}^{c_a}{p}_a^b)=\frac{{\Sigma }_{b=0}^{c_1+1}prime[j{]}^b}{{\Sigma }_{b=0}^{c_1}prime[j{]}^b}\ast {\Pi }_{a=1}^k({\Sigma }_{b=0}^{c_a}{p}_a^b)=\frac{{\Sigma }_{b=0}^{c_1+1}prime[j{]}^b}{{\Sigma }_{b=0}^{c_1}prime[j{]}^b}\ast sumd(i)$$
  

因此，要计算 $sumd(x)$ 的值，只需要记录其最小质因子的累计幂级数的和（${\Sigma }_{b=0}^{c_1}prime[j{]}^b$），依次递推，用与约数个数函数相同的方式计算即可。
代码如下。

#include <bits/stdc++.h> using namespace std; const int maxn = 100; long long sumd[maxn], pre[maxn]; int prime[maxn], cnt; bool vist[maxn + 1]; void Euler(int n){ 
    sumd[1] = 1; for (int i = 2; i <= n; i ++){ 
    if (!vist[i]) prime[++ cnt] = i, sumd[i] = 1 + i, pre[i] = 1 + i; for (int j = 1; j <= cnt && i * prime[j] <= n; j ++){ 
    vist[i * prime[j]] = true; if (i % prime[j] != 0){ 
    pre[i * prime[j]] = (1 + prime[j]); sumd[i * prime[j]] = sumd[i] * (1 + prime[j]); } else { 
    pre[i * prime[j]] = pre[i] * prime[j] + 1; sumd[i * prime[j]] = sumd[i] / pre[i] * pre[i * prime[j]]; break; } } } } int main(){ 
    Euler(maxn); for (int i = 1; i <= 100; i ++) printf("%d: %d\n", i, sumd[i]); return 0; }