一、问题引入
求 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(x⩾1) 倍都是合数,对于任意一个合数
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(x∈N,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])
φ(i∗prime[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) φ(i∗prime[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]
i∗prime[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) φ(i∗prime[j])=(i∗prime[j])∗Πi=1k(1−pi1)=prime[j]∗(i∗Πi=1k(1−pi1))=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(i∗prime[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]
i∗prime[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(i∗prime[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(i∗prime[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]
i∗prime[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(i∗prime[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]b∗sumd(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;
}