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

一、问题引入

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

二、问题分析

1. 朴素的素数筛法

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

2. Eratosthenes筛法

埃拉托斯特尼筛法，简称埃氏筛法、埃筛。
埃筛用到的原理十分简单：一个质数 $p$ 的任意 $x(x⩾1)$ 倍都是合数，对于任意一个合数 $c$ 都存在一个质数 $p$ 使得 $c=xp(x\in \mathbb{N},1<x<n)$。这句话的前半句保证了埃筛通过质数的若干倍筛去的数必定是合数，而后半句表示了埃筛必定可以在筛去所有的合数。采用从 $1$ 到 $n$ 的遍历顺序，就可以在遍历到某个合数之前筛去它，从而得到素数表。

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)$。实测运行时间大致如下。

数据范围	运行时间
$⩽10^6$	$⩽0.1s$
$⩽10^7$	$⩽0.2s$
$⩽10^8$	$\approx 1.5s$
需要说明的是，由于 $i^2$ 以下的合数必定被比 $i$ 小的质数筛掉了，因此每次筛选从 $i^2$ 开始到 $n$。

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);
}

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

数据范围	运行时间
$⩽10^6$	$⩽0.1s$
$⩽10^7$	$\approx 0.11s$
$⩽10^8$	$\approx 0.9s$

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筛法

即使对埃筛进行上述优化，依然会有合数会被多个质数筛掉。下面介绍Euler筛法，将筛法的复杂度优化至线性。
Euler筛法，又称欧拉筛、线性筛法、线筛。是一种可以在线性时间复杂度内求出 $1$ ~ $n$ 之间所有质数的筛法。思路是：让每个合数只被其最小的质因子筛去，从而保证每个合数只被筛去一次，做到线性复杂度。

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筛法求欧拉函数的值

关于欧拉函数的介绍，可以参考这一篇文章。
通过 $\phi (i)$ 和 $\phi (prime[j])$ 计算 $\phi (i\ast prime[j])$ 分成两种情况。

• 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)$$

又观察到：
$${\Sigma }_{b=0}^{c_1+1}prime[j{]}^b=({\Sigma }_{b=0}^{c_1}prime[j{]}^b)\ast prime[j]+1$$

因此，要计算 $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;
}