莫比乌斯反演是数论中的一种小技巧,常用于解决一些数论函数的求和问题。通常做法是使用 μ μ 函数来将 [n=1] [ n = 1 ] 转换成 ∑d|nμ(i) ∑ d | n μ ( i ) 的形式。
推荐博客:algocode 算法博客
μ μ 的定义:
- 如果 n n 有至少一个平方因子,那么
- 如果 n n 没有一个平方因子,且有 个素因子,那么 μ(n)=−1 μ ( n ) = − 1
- 如果 n n 没有一个平方因子,且有 个素因子,那么 μ(n)=1 μ ( n ) = 1
μ μ 的性质:
- μ μ 是积性函数。
- 根据这一性质,我们可以使用线性筛在 Θ(n) Θ ( n ) 的时间中筛出 i∈[1,n] i ∈ [ 1 , n ] 中的所有 μ(i) μ ( i ) 。
void prework() {
mu[1] = 1;
for (int i = 2; i <= 50000; i++) {
if (!mark[i]) {
prime[++tot] = i;
mu[i] = -1;
}
for (int j = 1; j <= tot && i * prime[j] <= 50000; j++) {
mark[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
}
mu[i * prime[j]] = -mu[i];
}
mu[i] += mu[i - 1];
}
}
- ∑d|nμ(i)=1 ∑ d | n μ ( i ) = 1
- 证明见 algocode 算法博客。
- 利用这一性质,我们就可以解决一些简单的问题了。
给定
a,b,c,d,e
a
,
b
,
c
,
d
,
e
,求
∑bi=a∑dj=c[gcd(i,j)=k]
∑
i
=
a
b
∑
j
=
c
d
[
g
c
d
(
i
,
j
)
=
k
]
的值。
1≤a≤b,c≤d,k≤50000
1
≤
a
≤
b
,
c
≤
d
,
k
≤
50000
,
50000
50000
组询问。
可以将询问拆成四个前缀询问。
问题变成了:求
∑ni=1∑nj=1[gcd(i,j)=k]
∑
i
=
1
n
∑
j
=
1
n
[
g
c
d
(
i
,
j
)
=
k
]
的值。
将 i,j i , j 都除以 k k ,上式等于
利用 μ μ 函数的第二个性质,上式等于 ∑⌊nk⌋i=1∑⌊nk⌋j=1∑d|i,d|jμ(d) ∑ i = 1 ⌊ n k ⌋ ∑ j = 1 ⌊ n k ⌋ ∑ d | i , d | j μ ( d )
将枚举
d
d
移到前面,并令 ,
则上式等于
∑nd=1μ(d)∑⌊nkd⌋i′=1∑⌊nkd⌋j′=11=∑nd=1μ(d)⌊nkd⌋2
∑
d
=
1
n
μ
(
d
)
∑
i
′
=
1
⌊
n
k
d
⌋
∑
j
′
=
1
⌊
n
k
d
⌋
1
=
∑
d
=
1
n
μ
(
d
)
⌊
n
k
d
⌋
2
而 ⌊ni⌋ ⌊ n i ⌋ 只有 Θ(n‾√) Θ ( n ) 种取值。直接预处理 μ μ 函数的前缀和,然后数论分块即可。
#include <cstdio>
#include <iostream>
using namespace std;
const int maxn = 50005;
bool mark[maxn];
int tot, prime[maxn], mu[maxn];
int T, a, b, c, d, k;
void prework() {
mu[1] = 1;
for (int i = 2; i <= 50000; i++) {
if (!mark[i]) {
prime[++tot] = i;
mu[i] = -1;
}
for (int j = 1; j <= tot && i * prime[j] <= 50000; j++) {
mark[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
}
mu[i * prime[j]] = -mu[i];
}
mu[i] += mu[i - 1];
}
}
int solve(int n, int m) {
int res = 0;
for (int i = 1, j; i <= min(n, m); i = j + 1) {
j = min(n / (n / i), m / (m / i));
res += (mu[j] - mu[i - 1]) * (n / i) * (m / i);
}
return res;
}
int main() {
prework();
scanf("%d", &T);
while (T--) {
scanf("%d %d %d %d %d", &a, &b, &c, &d, &k);
printf("%d\n", solve(b / k, d / k) - solve(b / k, (c - 1) / k) - solve((a - 1) / k, d / k) + solve((a - 1) / k, (c - 1) / k));
}
return 0;
}
求: ∑ni=1∑mj=1lcm(i,j) ∑ i = 1 n ∑ j = 1 m l c m ( i , j ) 。 n,m≤107 n , m ≤ 10 7
发现 lcm(i,j)=ijgcd(i,j) l c m ( i , j ) = i j g c d ( i , j ) 。
考虑枚举
gcd(i,j)=d
g
c
d
(
i
,
j
)
=
d
。问题就变成了:
令
f(n,m)
f
(
n
,
m
)
表示
∑ni=1∑mj=1[gcd(i,j)=1]⋅ij
∑
i
=
1
n
∑
j
=
1
m
[
g
c
d
(
i
,
j
)
=
1
]
⋅
i
j
,我们发现可以使用莫比乌斯反演得出原式等于:
这个式子可以使用数论分块来做。
然后,回到原问题。我们再使用一次数论分块,这样就可以在 Θ((n‾√)2)=Θ(n) Θ ( ( n ) 2 ) = Θ ( n ) 的时间内求解了。
#include <cstdio>
#include <iostream>
using namespace std;
typedef long long ll;
const int maxn = 10000005;
const int mod = 20101009;
bool mark[maxn];
ll n, m, tot, prime[maxn], mu[maxn], sum[maxn];
void prework() {
mu[1] = 1;
for (int i = 2; i <= min(n, m); i++) {
if (!mark[i]) {
prime[++tot] = i;
mu[i] = -1;
}
for (int j = 1; j <= tot && i * prime[j] <= min(n, m); j++) {
mark[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
}
mu[i * prime[j]] = -mu[i];
}
}
for (ll i = 1; i <= min(n, m); i++) {
sum[i] = (sum[i - 1] + i * i % mod * (mu[i] + mod)) % mod;
}
}
ll prod(ll x, ll y) {
return (x * (x + 1) / 2 % mod) * (y * (y + 1) / 2 % mod) % mod;
}
ll func(ll x, ll y) {
ll res = 0;
for (ll i = 1, j; i <= min(x, y); i = j + 1) {
j = min(x / (x / i), y / (y / i));
res = (res + (sum[j] - sum[i - 1] + mod) * prod(x / i, y / i)) % mod;
}
return res;
}
int main() {
scanf("%lld %lld", &n, &m);
prework();
ll ans = 0;
for (ll i = 1, j; i <= min(n, m); i = j + 1) {
j = min(n / (n / i), m / (m / i));
ans = (ans + (j - i + 1) * (i + j) / 2 % mod * func(n / i, m / i)) % mod;
}
printf("%lld\n", ans);
return 0;
}