筛法(埃筛、线性筛)学习笔记

本文介绍了求解闭区间内所有质数的两种经典方法:埃拉托斯特尼筛法和欧拉筛法,并详细分析了这两种筛法的实现原理、优化手段及时间复杂度。此外,还探讨了如何利用欧拉筛法求解特定数论函数的问题。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

一、问题引入

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

二、问题分析

1. 朴素的素数筛法

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

2. Eratosthenes筛法

埃拉托斯特尼筛法,简称埃氏筛法、埃筛。
埃筛用到的原理十分简单:一个质数 p p p 的任意 x ( x ⩾ 1 ) x(x \geqslant 1) x(x1) 倍都是合数,对于任意一个合数 c c c 都存在一个质数 p p p 使得 c = x p ( x ∈ N , 1 < x < n ) c=xp(x \in \mathbb{N}, 1\lt x\lt n) c=xp(xN,1<x<n)。这句话的前半句保证了埃筛通过质数的若干倍筛去的数必定是合数,而后半句表示了埃筛必定可以在筛去所有的合数。采用从 1 1 1 n n 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 ( n l o g l o g n ) O(nloglogn) O(nloglogn)。实测运行时间大致如下。

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

2.2 Eratosthenes筛法的时间效率优化

2.2.1 减小筛法上界

需要求出 1 1 1 n n n 的所有质数,并不需要取质数进行筛选,只需要取小于 n \sqrt{n} 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 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);
}

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

数据范围运行时间
⩽ 1 0 6 \leqslant 10^6 106 ⩽ 0.1 s \leqslant 0.1s 0.1s
⩽ 1 0 7 \leqslant 10^7 107 ≈ 0.11 s \approx 0.11s 0.11s
⩽ 1 0 8 \leqslant 10^8 108 ≈ 0.9 s \approx 0.9s 0.9s

2.3 Eratosthenes筛法的空间效率优化

2.3.1 质数个数的渐进

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

2.3.2 压位

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

2.3.3 分块存储

由于只需要取小于 n \sqrt{n} n 的质数进行筛选,因此,只需要保存 1 1 1 n \sqrt{n} n 的vist数组,然后对后面的数分成 ⌈ n s ⌉ \lceil \frac {n} {s}\rceil sn 个大小为 s s s 块,用前几个质数筛每一个块里的数,因此每次只需要 ⌈ n s ⌉ \lceil \frac {n} {s}\rceil sn 大小的数组,总空间复杂度为 O ( n + s ) O(\sqrt{n} + s) O(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 1 1 ~ n n 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 ) O(n) O(n),实测运行时间略慢于埃筛。这点令我十分惊讶,究其原因可能是因为Euler筛法需要多次取模,导致常数稍大。

3.2 Euler筛法的应用

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

3.2.1 Euler筛法求欧拉函数的值

关于欧拉函数的介绍,可以参考这一篇文章
通过 φ ( i ) \varphi(i) φ(i) φ ( p r i m e [ j ] ) \varphi(prime[j]) φ(prime[j]) 计算 φ ( i ∗ p r i m e [ j ] ) \varphi(i*prime[j]) φ(iprime[j]) 分成两种情况。

  • case 1: g c d ( i , p r i m e [ j ] ) = 1 gcd(i, prime[j]) = 1 gcd(i,prime[j])=1,由积性函数的性质有:
    φ ( i ∗ p r i m e [ j ] ) = φ ( i ) ∗ φ ( p r i m e [ j ] ) = φ ( i ) ∗ ( p r i m e [ j ] − 1 ) \varphi(i*prime[j]) = \varphi(i) * \varphi(prime[j]) = \varphi(i) * (prime[j] - 1) φ(iprime[j])=φ(i)φ(prime[j])=φ(i)(prime[j]1)
  • case 2: g c d ( i , p r i m e [ j ] ) = p r i m e [ j ] gcd(i, prime[j]) = prime[j] gcd(i,prime[j])=prime[j],根据 i ∗ p r i m e [ j ] i*prime[j] iprime[j] i i i 的唯一分解、欧拉函数的性质有:
    φ ( i ∗ p r i m e [ j ] ) = ( i ∗ p r i m e [ j ] ) ∗ Π i = 1 k ( 1 − 1 p i ) = p r i m e [ j ] ∗ ( i ∗ Π i = 1 k ( 1 − 1 p i ) ) = p r i m e [ j ] ∗ φ ( i ) \varphi(i*prime[j]) = (i*prime[j])*\Pi_{i=1}^{k}(1-\frac{1}{p_i})=prime[j]*(i*\Pi_{i=1}^{k}(1-\frac{1}{p_i})) = prime[j] * \varphi(i) φ(iprime[j])=(iprime[j])Πi=1k(1pi1)=prime[j](iΠi=1k(1pi1))=prime[j]φ(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 ) d(x) d(x) 表示,其值为 x x x 的约数个数。根据 x x x 的唯一分解和乘法原理可得:
x = Π a = 1 k p a c a x = \Pi_{a = 1}^{k}p_a^{c_a} x=Πa=1kpaca,则 d ( x ) = Π i = a k ( c a + 1 ) d(x) = \Pi_{i =a}^{k}(c_a + 1) d(x)=Πi=ak(ca+1)
易证约数个数函数为积性函数,尝试通过Euler筛法计算,分类讨论。

  • case 1: g c d ( i , p r i m e [ j ] ) = 1 gcd(i, prime[j]) = 1 gcd(i,prime[j])=1,由积性函数性质有:
    d ( i ∗ p r i m e [ j ] ) = d ( i ) ∗ d ( p r i m e [ j ] ) = 2 d ( i ) d(i * prime[j]) = d(i) * d(prime[j]) = 2d(i) d(iprime[j])=d(i)d(prime[j])=2d(i)
  • case 2: g c d ( i , p r i m e [ j ] ) = p r i m e [ j ] gcd(i, prime[j]) = prime[j] gcd(i,prime[j])=prime[j],由 i ∗ p r i m e [ j ] i*prime[j] iprime[j] i i i 的唯一分解、约数个数函数的性质有(由于 p r i m e [ j ] prime[j] prime[j] i i i 的最小质因子,不妨设 i = ( p r i m e [ j ] ) c 1 ∗ Π a = 2 k p a c a i = (prime[j])^{c_1}*\Pi_{a=2}^{k}p_a^{c_a} i=(prime[j])c1Πa=2kpaca):
    d ( i ∗ p r i m e [ j ] ) = ( c 1 + 2 ) ∗ Π a = 2 k ( c a + 1 ) = c 1 + 2 c 1 + 1 ∗ Π a = 1 k c a = c 1 + 2 c 1 + 1 d ( i ) d(i * prime[j]) = (c_1+2)*\Pi_{a=2}^{k}(c_a + 1) = \frac{c_1+2}{c_1+1}*\Pi_{a=1}^{k}c_a = \frac{c_1+2}{c_1+1}d(i) d(iprime[j])=(c1+2)Πa=2k(ca+1)=c1+1c1+2Πa=1kca=c1+1c1+2d(i)

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

#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筛法求约数和函数的值

约束和函数,记作 s u m d ( x ) sumd(x) sumd(x),其值为 x x x 的所有约数的和。根据 x x x 的唯一分解和加乘原理可得:
x = Π a = 1 k p a c a x = \Pi_{a=1}^{k}p_a^{c_a} x=Πa=1kpaca,则 s u m d ( x ) = Π a = 1 k ( Σ b = 0 c a p a b ) sumd(x) = \Pi_{a=1}^{k}(\Sigma_{b=0}^{c_a}p_a^{b}) sumd(x)=Πa=1k(Σb=0capab)
同样易证,约束和函数为积性函数,尝试通过Euler筛法计算,分类讨论。

  • case 1: g c d ( i , p r i m e [ j ] ) = 1 gcd(i, prime[j])=1 gcd(i,prime[j])=1,由积性函数性质有:
    s u m d ( i ∗ p r i m e [ j ] ) = s u m d ( i ) ∗ s u m d ( p r i m e [ j ] ) = s u m ( d ) ∗ ( 1 + p r i m e [ j ] ) sumd(i * prime[j]) = sumd(i) * sumd(prime[j]) = sum(d) * (1 + prime[j]) sumd(iprime[j])=sumd(i)sumd(prime[j])=sum(d)(1+prime[j])
  • case 2: g c d ( i , p r i m e [ j ] ) = p r i m e [ j ] gcd(i, prime[j]) = prime[j] gcd(i,prime[j])=prime[j],由 i ∗ p r i m e [ j ] i*prime[j] iprime[j] i i i 的唯一分解、约数个数函数的性质有(由于 p r i m e [ j ] prime[j] prime[j] i i i 的最小质因子,不妨设 i = ( p r i m e [ j ] ) c 1 ∗ Π a = 2 k p a c a i = (prime[j])^{c_1}*\Pi_{a=2}^{k}p_a^{c_a} i=(prime[j])c1Πa=2kpaca):
    s u m d ( i ∗ p r i m e [ j ] ) = Σ b = 0 c 1 + 1 p r i m e [ j ] b ∗ Π a = 2 k ( Σ b = 0 c a p a b ) = Σ b = 0 c 1 + 1 p r i m e [ j ] b Σ b = 0 c 1 p r i m e [ j ] b ∗ Π a = 1 k ( Σ b = 0 c a p a b ) = Σ b = 0 c 1 + 1 p r i m e [ j ] b Σ b = 0 c 1 p r i m e [ j ] b ∗ s u m d ( i ) sumd(i * prime[j]) = \Sigma_{b=0}^{c_1 + 1}prime[j]^{b}*\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}}*\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}}*sumd(i) sumd(iprime[j])=Σb=0c1+1prime[j]bΠa=2k(Σb=0capab)=Σb=0c1prime[j]bΣb=0c1+1prime[j]bΠa=1k(Σb=0capab)=Σb=0c1prime[j]bΣb=0c1+1prime[j]bsumd(i)

又观察到:
Σ b = 0 c 1 + 1 p r i m e [ j ] b = ( Σ b = 0 c 1 p r i m e [ j ] b ) ∗ p r i m e [ j ] + 1 \Sigma_{b=0}^{c_1 + 1}prime[j]^{b} = (\Sigma_{b=0}^{c_1 }prime[j]^{b}) * prime[j] + 1 Σb=0c1+1prime[j]b=(Σb=0c1prime[j]b)prime[j]+1

因此,要计算 s u m d ( x ) sumd(x) sumd(x) 的值,只需要记录其最小质因子的累计幂级数的和( Σ b = 0 c 1 p r i m e [ j ] b \Sigma_{b=0}^{c_1 }prime[j]^{b} Σb=0c1prime[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;
}
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值