本文详细介绍了线性筛和埃氏筛两种算法，用于高效地计算质数、欧拉函数和莫比乌斯函数。线性筛能在线性时间内筛出质数，而埃氏筛则常用于计算欧拉函数和莫比乌斯函数。通过实例代码展示了这两种筛法的具体实现，并给出了在实际问题中如P2185仪仗队和P2522[HAOI2011]Problem B中的应用，展示了如何利用筛法解决数论问题。
传送门：数论函数学习笔记
埃氏筛
埃氏筛的思想很简单，就是从前往后，遇到一个质数就把所有后面小于 n 的这个质数的倍数全都标记为合数。后面遇到他们就直接跳过。
一、筛质数
按照上面的思想直接模拟就好了，时间复杂度 O ( n log log n ) O(n \log \log n) O(nloglogn)：
#include<bits/stdc++.h>
using namespace std;
#define MAXN 100100
int v[MAXN] = { 0 }; // 合数标记
void primes(int n){
memset(v, 0, sizeof(v));
for(int i = 2; i <= n; i++){
if(v[i]) continue;
cout << i << endl;
for(int j = 1; j <= n / i; j++) v[i * j] = 1;
}
}
int main(){
primes(100);
return 0;
}
二、筛欧拉函数
利用欧拉函数的计算式和埃氏筛，能在 O ( n log n ) O(n \log n) O(nlogn) 的时间内筛出欧拉函数的值。就是首先初始化所有的 phi[i] = 1，然后再枚举从 2 到 n 的每个数，用埃氏筛找出质数，然后把后面质数的倍数根据计算式计算就好了。
#include<bits/stdc++.h>
using namespace std;
#define MAXN 100100
int phi[MAXN] = { 0 };
void euler(int n){
for(int i = 1; i <= n; i++) phi[i] = i;
for(int i = 2; i <= n; i++){
if(phi[i] == i)
for(int j = i; j <= n; j += i)
phi[j] = phi[j] / i * (i - 1);
}
}
int main(){
euler(20);
for(int i = 1; i <= 20; i++){
cout << "phi[" << i << "] = " << phi[i] << endl;
}
return 0;
}
三、筛莫比乌斯函数
首先把所有的 mu 初始化为 1，接下来对于每一个质数 p 把 mu[p] 变成 -1。并扫描 p 的倍数们看他们能否被 p 2 p^2 p2 整除，如果可以 mu[x] = 0，如果不行 mu[x] = -mu[x]。
#include<bits/stdc++.h>
using namespace std;
#define MAXN 100100
int v[MAXN] = { 0 };
int mu[MAXN] = { 0 };
void mobius(int n){
memset(v, 0, sizeof(v));
for(int i = 2; i <= n; i++) mu[i] = 1;
for(int i = 2; i <= n; i++){
if(v[i]) continue;
mu[i] = -1;
for(int j = 2 * i; j <= n; j += i){
v[j] = 1;
if((j / i) % i == 0) mu[j] = 0;
else mu[j] *= -1;
}
}
}
int main(){
mobius(20);
for(int i = 1; i <= 20; i++){
cout << "mu[" << i << "] = " << mu[i] << endl;
}
return 0;
}
线性筛
一、筛质数
顾名思义，这玩意儿能在线性的时间内筛出质数。具体的做法是优化一下埃氏筛，在埃氏筛里面，一个合数可能被很多个质数标记过，比如 12 就被 2 和 3 同时标记过。这样就会降低我们的效率。线性筛的想法就是让所有的合数只有一种被筛出来的方式。也就是被它的最小质因子筛出来。具体做法如下：
- 枚举 2 ∼ n 2 \sim n 2∼n 的所有数，考虑用一下方式维护一个数组 v（表示这个数的最小质因子）。
- 如果 v[i] = i，那么这个是是质数，把它存下来。
- 扫描不大于 v[i] 的所有质数 p，令 v[ip] = p（很好理解吧）。
#include<bits/stdc++.h>
using namespace std;
#define MAXN 100100
int v[MAXN] = { 0 };
int prime[MAXN] = { 0 };
int primes(int n){
memset(v, 0, sizeof(v));
int cnt = 0; // 质数的个数
for(int i = 2; i <= n; i++){
if(!v[i]){
v[i] = i; prime[++cnt] = i;
}
for(int j = 1; j <= cnt; j++){ // 枚举所有质数
if(prime[j] > v[i] or prime[j] > n / i) break; // 质数比它大或者超过范围
v[i * prime[j]] = prime[j]; // 更新最小质因子
}
}
return cnt;
}
int main(){
int num = primes(100);
for(int i = 1; i <= num; i++){
cout << prime[i] << endl;
}
return 0;
}
二、筛欧拉函数
这就要用到几个欧拉函数的性质了：
- 如果 p ∣ n p \mid n p∣n 且 p 2 ∣ n p^2 \mid n p2∣n 则有 φ ( n ) = φ ( n p ) × p \varphi(n) = {\varphi(\frac np)} \times p φ(n)=φ(pn)×p
- 如果 p ∣ n p \mid n p∣n 且 p 2 ∤ n p^2 \nmid n p2∤n 则有 φ ( n ) = φ ( n p ) × ( p − 1 ) \varphi(n) = \varphi(\frac np) \times (p-1) φ(n)=φ(pn)×(p−1)
我们发现这两个判断和线性筛的两个判断很类似，所以我们可以利用线性筛，从 φ ( n p ) \varphi(\frac np) φ(pn) 递推到 φ ( n ) \varphi(n) φ(n)
#include<bits/stdc++.h>
using namespace std;
#define MAXN 100100
int v[MAXN] = { 0 };
int phi[MAXN] = { 0 };
void euler(int n){
memset(v, 0, sizeof(v));
int cnt = 0;
for(int i = 2; i <= n; i++){
if(!v[i]){
v[i] = i; prime[++cnt] = i;
phi[i] = i - 1;
}
for(int j = 1; j <= cnt; j++){
if(prime[j] > v[i] or prime[j] > n / i) break;
v[i *prime[j]] = prime[j];
if(i % prime[j]) phi[i * prime[j]] = phi[i] * (prime[j] - 1);
else phi[i * prime[j]] = phi[i] * prime[j];
}
}
}
int main(){
euler(20);
for(int i = 1; i <= 20; i++){
cout << "phi[" << i << "] = " << phi[i] << endl;
}
return 0;
}
P2185 仪仗队
题目大意是这样的，有一个点阵（方阵），你站在左下角，问你能总共看到几个点。然后我们可以发现，如果我们将左下角的点放在坐标轴 原点，能被看见的点的横纵坐标一定是互质的（除了两个和它相邻的点）。所以我们要求解的问题的答案就是：
a
n
s
=
2
+
∑
i
=
1
n
∑
j
=
1
n
[
g
c
d
(
i
,
j
)
=
1
]
ans = 2 + \sum_{i=1}^n \sum_{j=1}^n [gcd(i, j) = 1]
ans=2+i=1∑nj=1∑n[gcd(i,j)=1]
然后可以对这个式子进行化简：
a
n
s
−
2
=
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
[
g
c
d
(
i
,
j
)
=
1
]
=
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
ε
(
g
c
d
(
i
,
j
)
)
=
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
∑
d
∣
g
c
d
(
i
,
j
)
μ
(
d
)
=
∑
d
=
1
n
−
1
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
[
d
∣
g
c
d
(
i
,
j
)
]
μ
(
d
)
=
∑
d
=
1
n
−
1
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
[
d
∣
i
∧
d
∣
j
]
μ
(
d
)
=
∑
d
=
1
n
−
1
μ
(
d
)
∑
i
=
1
n
−
1
∑
j
=
1
n
−
1
[
d
∣
i
∧
d
∣
j
]
=
∑
d
=
1
n
−
1
μ
(
d
)
⌊
n
−
1
d
⌋
2
\begin{aligned} ans-2 = & \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} [gcd(i, j) = 1] \\ = & \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} \varepsilon(gcd(i, j)) = \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} \sum_{d \mid gcd(i, j)}\mu(d) \\ = & \sum_{d = 1}^{n-1} \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} [d \mid gcd(i, j)]\mu(d) = \sum_{d = 1}^{n-1} \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} [d \mid i \wedge d \mid j]\mu(d) \\ = & \sum_{d=1}^{n-1} \mu(d) \sum_{i=1}^{n-1} \sum_{j=1}^{n-1} [d \mid i \wedge d \mid j] \\ = & \sum_{d=1}^{n-1} \mu(d) \left\lfloor \frac {n-1}d \right\rfloor^2 \end{aligned}
ans−2=====i=1∑n−1j=1∑n−1[gcd(i,j)=1]i=1∑n−1j=1∑n−1ε(gcd(i,j))=i=1∑n−1j=1∑n−1d∣gcd(i,j)∑μ(d)d=1∑n−1i=1∑n−1j=1∑n−1[d∣gcd(i,j)]μ(d)=d=1∑n−1i=1∑n−1j=1∑n−1[d∣i∧d∣j]μ(d)d=1∑n−1μ(d)i=1∑n−1j=1∑n−1[d∣i∧d∣j]d=1∑n−1μ(d)⌊dn−1⌋2
然后 μ ( d ) \mu(d) μ(d) 可以用埃氏筛预处理然后 O ( 1 ) O(1) O(1) 查询。后面的 ⌊ n d ⌋ 2 \left\lfloor \frac nd \right\rfloor^2 ⌊dn⌋2 显然也是可以 O ( 1 ) O(1) O(1) 查询的，那么总的复杂度就是 O ( n log n + n ) O(n\log n + n) O(nlogn+n)
#include<bits/stdc++.h>
using namespace std;
#define MAXN 40400
#define endl '\n'
int n = 0;
int v[MAXN] = { 0 };
int mu[MAXN] = { 0 };
void mobius(int n){
memset(v, 0, sizeof(v));
for(int i = 1; i <= n; i++) mu[i] = 1;
for(int i = 2; i <= n; i++){
if(v[i]) continue;
mu[i] = -1;
for(int j = 2 * i; j <= n; j += i){
v[j] = 1;
if((j / i) % i == 0) mu[j] = 0;
else mu[j] *= -1;
}
}
}
int main(){
scanf("%d", &n); n -= 1;
if(n == 0){
cout << 0 << endl; return 0;
}
mobius(n);
int ans = 0;
for(int d = 1; d <= n; d++){
int temp = mu[d] * (n / d) * (n / d);
ans += temp;
}
ans += 2;
cout << ans << endl;
return 0;
}
P2522 [HAOI2011]Problem b
P2522 [HAOI2011]Problem b
题目大意：对于给出的 n 个询问，每次求出有多少个数对
(
x
,
y
)
(x, y)
(x,y) 满足
a
≤
x
≤
b
,
c
≤
y
≤
d
a \leq x \leq b, c \leq y \leq d
a≤x≤b,c≤y≤d 且
g
c
d
(
x
,
y
)
=
k
gcd(x, y) = k
gcd(x,y)=k。
所以这个题的答案就是：
a
n
s
=
∑
i
=
a
b
∑
j
=
c
d
[
g
c
d
(
i
,
j
)
=
k
]
ans = \sum_{i = a}^b \sum_{j = c}^d [gcd(i, j) = k]
ans=i=a∑bj=c∑d[gcd(i,j)=k]
这个式子就和上面那个题的试着长得很像了，但是这个上下界看着还是有点不舒服，因为我们想把它转化成下界为 1 的求和，于是我们利用容斥原理，转化一下问题：
a
n
s
=
∑
i
=
a
b
∑
j
=
c
d
[
g
c
d
(
i
,
j
)
=
k
]
=
∑
i
=
1
b
∑
j
=
1
d
[
g
c
d
(
i
,
j
)
=
k
]
−
∑
i
=
1
a
−
1
∑
j
=
1
d
[
g
c
d
(
i
,
j
)
=
k
]
−
∑
i
=
1
b
∑
j
=
1
c
−
1
[
g
c
d
(
i
,
j
)
=
k
]
+
∑
i
=
1
a
−
1
∑
j
=
1
c
−
1
[
g
c
d
(
i
,
j
)
=
k
]
\begin{aligned} ans = & \sum_{i = a}^b \sum_{j = c}^d [gcd(i, j) = k] \\ = & \sum_{i = 1}^b \sum_{j = 1}^d [gcd(i, j) = k] - \sum_{i = 1}^{a-1} \sum_{j = 1}^d [gcd(i, j) = k] - \sum_{i = 1}^b \sum_{j = 1}^{c-1} [gcd(i, j) = k] + \sum_{i = 1}^{a-1} \sum_{j = 1}^{c-1} [gcd(i, j) = k] \end{aligned}
ans==i=a∑bj=c∑d[gcd(i,j)=k]i=1∑bj=1∑d[gcd(i,j)=k]−i=1∑a−1j=1∑d[gcd(i,j)=k]−i=1∑bj=1∑c−1[gcd(i,j)=k]+i=1∑a−1j=1∑c−1[gcd(i,j)=k]
于是问题就被转化成求 4 个和上道题差不多的式子。由上一道题我们知道：
∑
i
=
1
n
∑
j
=
1
m
[
g
c
d
(
i
,
j
)
=
k
]
=
∑
i
=
1
n
∑
j
=
1
m
ε
(
g
c
d
(
i
,
j
)
k
)
=
∑
i
=
1
n
∑
j
=
1
m
∑
p
∣
g
c
d
(
i
,
j
)
/
k
μ
(
p
)
=
∑
p
=
1
m
i
n
(
n
,
m
)
∑
i
=
1
n
∑
j
=
1
m
[
p
∣
g
c
d
(
i
,
j
)
k
]
μ
(
p
)
=
∑
p
=
1
m
i
n
(
n
,
m
)
μ
(
p
)
∑
i
=
1
n
∑
j
=
1
m
[
p
∣
g
c
d
(
i
,
j
)
k
]
=
∑
p
=
1
m
i
n
(
n
,
m
)
μ
(
p
)
∑
i
=
1
n
∑
j
=
1
m
[
p
∣
i
k
∧
p
∣
j
k
]
=
∑
p
=
1
m
i
n
(
n
,
m
)
μ
(
p
)
∑
i
=
1
n
∑
j
=
1
m
[
k
p
∣
i
∧
k
p
∣
j
]
=
∑
p
=
1
m
i
n
(
n
,
m
)
μ
(
p
)
⌊
n
k
p
⌋
⌊
m
k
p
⌋
\begin{aligned} & \sum_{i = 1}^n \sum_{j = 1}^m [gcd(i, j) = k]\\ = & \sum_{i = 1}^n \sum_{j = 1}^m \varepsilon\left( \frac{gcd(i, j)}{k} \right)= \sum_{i = 1}^n \sum_{j = 1}^m \sum_{p | gcd(i, j) / k} \mu(p) = \sum_{p = 1}^{min(n, m)}\sum_{i = 1}^n \sum_{j = 1}^m \left[ p \Big| \frac{gcd(i, j)}{k} \right]\mu(p) \\ = &\sum_{p = 1}^{min(n, m)}\mu(p) \sum_{i = 1}^n \sum_{j = 1}^m \left[ p \Big| \frac{gcd(i, j)}{k} \right] = \sum_{p = 1}^{min(n, m)}\mu(p) \sum_{i = 1}^n \sum_{j = 1}^m \left[ p \Big| \frac{i}{k} \wedge p \Big| \frac{j}{k} \right] = \sum_{p = 1}^{min(n, m)}\mu(p) \sum_{i = 1}^n \sum_{j = 1}^m \left[ kp \mid i \wedge kp \mid j \right] \\ = & \sum_{p = 1}^{min(n, m)}\mu(p) \left\lfloor \frac{n}{kp} \right\rfloor \left\lfloor \frac{m}{kp} \right\rfloor \end{aligned}
===i=1∑nj=1∑m[gcd(i,j)=k]i=1∑nj=1∑mε(kgcd(i,j))=i=1∑nj=1∑mp∣gcd(i,j)/k∑μ(p)=p=1∑min(n,m)i=1∑nj=1∑m[p∣∣∣kgcd(i,j)]μ(p)p=1∑min(n,m)μ(p)i=1∑nj=1∑m[p∣∣∣kgcd(i,j)]=p=1∑min(n,m)μ(p)i=1∑nj=1∑m[p∣∣∣ki∧p∣∣∣kj]=p=1∑min(n,m)μ(p)i=1∑nj=1∑m[kp∣i∧kp∣j]p=1∑min(n,m)μ(p)⌊kpn⌋⌊kpm⌋
然后我们只需要写一个函数来处理这个，然后每次询问调用 4 次就好了。计算这个式子的时间复杂度是 O ( m i n ( n , m ) log m i n ( n , m ) + m i n ( n , m ) ) O(min(n, m) \log min(n, m) + min(n, m)) O(min(n,m)logmin(n,m)+min(n,m))，所以如果我还要处理 5 × 1 0 4 5 \times 10^4 5×104 组这样的数据，就会超时，直接计算能拿 30 分。
#include<bits/stdc++.h>
using namespace std;
#define MAXN 50500
#define in read()
inline int read(){
int x = 0; char c = getchar();
while(c < '0' or c > '9') c = getchar();
while('0' <= c and c <= '9'){
x = x * 10 + c - '0'; c = getchar();
}
return x;
}
int t = 0; int k = 0;
int a = 0; int b = 0;
int c = 0; int d = 0;
int v[MAXN] = { 0 };
int mu[MAXN] = { 0 };
void mobius(int n){
memset(v, 0, sizeof(v));
for(int i = 1; i <= n; i++) mu[i] = 1;
for(int i = 2; i <= n; i++){
if(v[i]) continue;
mu[i] = -1;
for(int j = 2 * i; j <= n; j += i){
v[j] = 1;
if((j / i) % i == 0) mu[j] = 0;
else mu[j] *= -1;
}
}
}
int cal(int n, int m){
int r = min(n, m);
int ans = 0;
for(int i = 1; i <= r; i++){
int temp = mu[i] * (n / (i * k)) * (m / (i * k));
ans += temp;
}
return ans;
}
int main(){
t = in;
mobius(50000);
while(t--){
a = in; b = in;
c = in; d = in;
k = in;
int ans = cal(b, d) - cal(a - 1, d) - cal(b, c - 1) + cal(a - 1, c - 1);
cout << ans << '\n';
}
return 0;
}
下面我们考虑怎么样优化处理这个式子的时间。首先看看这个式子的样子，他长得是不是就很像整除分块qwq，那就直接整除分块计算这个式子，复杂度就被降到了 O ( m i n ( n , m ) ) O(\sqrt{min(n, m)}) O(min(n,m))，·这样一来就能过了（不知道为什在 luogu 上过不了，但是在 bzoj上就一次过了）。bzoj 传送门：bzoj 2301
#include<bits/stdc++.h>
using namespace std;
#define MAXN 50500
#define in read()
inline int read(){
int x = 0; char c = getchar();
while(c < '0' or c > '9') c = getchar();
while('0' <= c and c <= '9'){
x = x * 10 + c - '0'; c = getchar();
}
return x;
}
int t = 0; int k = 0;
int a = 0; int b = 0;
int c = 0; int d = 0;
int v[MAXN] = { 0 };
int mu[MAXN] = { 0 };
int sum[MAXN] = { 0 };
void mobius(int n){
for(int i = 1; i <= n; i++) mu[i] = 1;
for(int i = 2; i <= n; i++){
if(v[i]) continue;
mu[i] = -1;
for(int j = 2 * i; j <= n; j += i){
v[j] = 1;
if((j / i) % i == 0) mu[j] = 0;
else mu[j] *= -1;
}
}
for(int i = 1; i <= n; i++) sum[i] = sum[i - 1] + mu[i];
}
int cal(int n, int m){
int x = min(n, m);
int ans = 0;
for(int l = 1, r; l <= x; l = r + 1){
r = min(n / (n / l), m / (m / l));
ans += (n / (l * k)) * (m / (l * k)) * (sum[r] - sum[l - 1]);
}
return ans;
}
int main(){
t = in;
mobius(50000);
while(t--){
a = in; b = in;
c = in; d = in;
k = in;
int ans = cal(b, d) - cal(a - 1, d) - cal(b, c - 1) + cal(a - 1, c - 1);
cout << ans << '\n';
}
return 0;
}
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
