莫比乌斯反演
/ / 阅读耗时 27 分钟 莫比乌斯反演及常见套路，不涉及证明。莫比乌斯反演难在公式推导和应用。
整除分块
在介绍莫比乌斯反演之前，先了解整除分块。
整除分块需要解决的问题是，如何求解下面的求和式：
当然可以$O(n)$解决，但那样在某些问题中不是正确复杂度。注意求和元素有许多重复的部分，它们呈块状分布，并且可以发现每一个块中，最右边的值为n/(n/l)。
那么就可以得到下面的快速计算算法：1
2
3
4
5int solve(int x) {
int l = 1, r, ans = 0;
while (l <= x)r = x / (x / l), ans += (r - l + 1) * (x / l), l = r + 1;
return ans;
}
这就是整除分块算法，可以在$O(\sqrt n)$下完成上面求和式的计算。有时候右边的求和元素可能乘一个函数，这时维护一下前缀和即可。
莫比乌斯函数
莫比乌斯函数$\mu(x)$是莫比乌斯反演中的重要函数，它是定义在正整数集下的函数，具体定义如下：
- $\mu(1)=1$
- 当$x=\displaystyle\prod_{i=1}^s p_i^{\alpha_i}$时（这里p为质数），如果所有的质数指数$\alpha_i$均为1，则$\mu(x)=(-1)^s$，否则为0。
莫比乌斯函数可以通过线性筛筛出来，类似筛欧拉函数的过程，只需将欧拉筛修改一下即可：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16int mu[N+1] = {0};//记录莫比乌斯函数值
int mark[N+1] = {0};//记录是不是素数
int prim[N+1], tot = 0;//储存当前已知的素数,tot记录个数
for (int i = 2; i <= N; i++) {
if (!mark[i]) {
prim[++tot] = i;
mu[i] = -1;//判断为素数,直接给莫比乌斯函数赋值
}
for (int j = 1; j <= tot; j++) {
if (i * prim[j] > N)break;
mark[i * prim[j]] = 1;//标记这个数一定不是素数
if (i % prim[j] == 0)break;//其实应该赋0，但是初始化即为0，省略这一步
else mu[i * prim[j]] = -mu[i];
}
}
别忘了mu[1]=1。
莫比乌斯函数有以下性质：
- 对于任何正整数n，有$\displaystyle\sum_{d|n}\mu(d)=[n=1]$，这里[P]在P条件成立时为1，否则为0。这个性质很常用。
- 莫比乌斯函数与欧拉函数的关系：$\displaystyle\sum_{d|n} \mu(d)\frac{n} {d}=\phi(n)$。
- 莫比乌斯函数是积性函数，即在$gcd(a,b)=1$时，有$\mu(ab)=\mu(a)\mu(b)$。
莫比乌斯反演定理
对于两个定义在自然数集上的函数$F(x)$和$f(x)$，若有：
那么有：
这个定理称为莫比乌斯反演定理。
莫比乌斯反演定理还有另一种形式，当$F(x)$和$G(x)$满足：
那么有：
关于公式推导
莫比乌斯反演的应用难在推导，这里介绍几个必要的公式。
这是一个比较显然的结论。
这个式子也有简化运算的功能，证明如下：
假设$x=kab+p,0\leq p<ab$，则$\lfloor \frac {x} {ab} \rfloor=k$。
那么可以推得$\lfloor \frac {x} {a}\rfloor=kb+\lfloor \frac {p} {a} \rfloor$。从而$\lfloor \frac {\frac {x} {a}} {b}\rfloor=k+\lfloor \frac {\lfloor \frac {p} {a} \rfloor} {b}\rfloor$。由于$p<ab$，故后一项必然为0。得证。
这个式子有一个推论：如果$\lfloor \frac {n} {a}\rfloor=\lfloor \frac {n} {b}\rfloor$，那么有$\lfloor \frac {n} {pa}\rfloor=\lfloor \frac {n} {pb}\rfloor$，这里p是正整数。
其余以后更新。
经典例题
介绍几个例题，不涉及题面，可点开标题链接自查。
ZAP-Queries
等价于求：
$[P]$在P条件成立时为1，否则为0。
关于gcd的莫比乌斯反演题目，大多走两个套路，这里是第一种：设$\displaystyle f(x)=\sum_{i=1}^a\sum_{j=1}^b[gcd(i,j)=x]$，并设$\displaystyle g(x)=\sum_{i=1}^a\sum_{j=1}^b[gcd(i,j)=sx]$，s为正整数。容易发现，有下面的关系：
莫比乌斯反演一步得到：
注意到g(x)是可以直接求出来的：
那么f(x)就是：
答案即为：
枚举$\frac {y} {d}$，可以得到：
这里需要满足$a\leq b$。到这一步就可以计算了，后面的求和式用整除分块去计算。时间复杂度$O(\sqrt n)$，相当高效。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
using namespace std;
typedef long long ll;
int mu[N + 1], prim[N + 1], mark[N + 1], tot, n, m, k;
inline int read() {
char e = getchar();
int s = 0;
while (e < '-')e = getchar();
while (e > '-')s = (s << 1) + (s << 3) + (e & 15), e = getchar();
return s;
}
inline void getMu() {//求莫比乌斯函数
for (int i = 2; i <= N; i++) {
if (!mark[i])prim[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot; j++) {
if ((ll) i * prim[j] > N)break;
mark[i * prim[j]] = 1;
if (i % prim[j] == 0)break;
else mu[i * prim[j]] = -mu[i];
}
}
mu[1] = 1;
}
int main() {
int t = read();
getMu();
for (int i = 1; i <= N; i++)mu[i] += mu[i - 1];//维护前缀和
while (t--) {
n = read(), m = read(), k = read();
if (n > m)swap(n, m);
int l = 1, r;
ll ans = 0;
while (l <= n) {
r = min(n / (n / l), m / (m / l));
ans += (ll) (n / (l * k)) * (m / (l * k)) * (mu[r] - mu[l - 1]);
l = r + 1;
}
printf("%lld\n", ans);
}
return 0;
}
YY的GCD
其实就是求：
直接用上一个问题的结论，得到：
枚举ip，得到：
就是：
后一个求和式可以用埃氏筛法的形式筛出来，然后整除分块处理。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
using namespace std;
typedef long long ll;
int mu[N + 1], prim[N + 1], mark[N + 1], tot, n, m;
ll sum[N + 1];
inline int read() {
char e = getchar();
int s = 0;
while (e < '-')e = getchar();
while (e > '-')s = (s << 1) + (s << 3) + (e & 15), e = getchar();
return s;
}
inline void getMu() {
for (int i = 2; i <= N; i++) {
if (!mark[i])prim[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot; j++) {
if ((ll) i * prim[j] > N)break;
mark[i * prim[j]] = 1;
if (i % prim[j] == 0)break;
else mu[i * prim[j]] = -mu[i];
}
}
mu[1] = 1;
}
int main() {
int t = read();
getMu();
for (int i = 1; i <= tot; i++) {
for (int j = 1; (ll) prim[i] * j < N; j++)sum[prim[i] * j] += mu[j];//筛后一个求和式的过程
}
for (int i = 1; i < N; i++)sum[i] += sum[i - 1];
while (t--) {
n = read(), m = read();
if (n > m)swap(n, m);
int l = 1, r;
ll ans = 0;
while (l <= n) {
r = min(n / (n / l), m / (m / l));
ans += (ll) (n / l) * (m / l) * (sum[r] - sum[l - 1]);
l = r + 1;
}
printf("%lld\n", ans);
}
return 0;
}
约数个数和
已经描述的很明确了，求：
$d(x)$为x的约数个数，这里有一个公式：
其实我并不会证，然后原式即为：
gcd问题第二个套路出现了：利用莫比乌斯函数性质得到：
枚举d：
提出d所在的求和：
枚举因子，这个操作很常用：
只枚举i、j关于d的倍数，得到：
就是：
这样仍然不容易看出算法来，再来一步：
后面两个求和都是整除分块的前缀和，预处理一下，再用整除分块求和即可。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
using namespace std;
typedef long long ll;
inline int read() {
char e = getchar();
int s = 0;
while (e < '-')e = getchar();
while (e > '-')s = (s << 1) + (s << 3) + (e & 15), e = getchar();
return s;
}
ll s[N + 5];
int mu[N + 5], prim[N + 5], tot, mark[N + 5];
inline void init() {
for (int i = 2; i <= N; i++) {
if (!mark[i]) {
prim[++tot] = i;
mu[i] = -1;
}
for (int j = 1; j <= tot; j++) {
if (prim[j] * (ll) i > N)break;
mark[i * prim[j]] = 1;
if (i % prim[j] == 0)break;
mu[i * prim[j]] = -mu[i];
}
}
mu[1] = 1;
for (int i = 2; i <= N; i++)mu[i] += mu[i - 1];//前缀和
for (int i = 1; i <= N; i++) {//预处理
int l = 1, r;
while (l <= i)r = i / (i / l), s[i] += (ll) (r - l + 1) * (i / l), l = r + 1;
}
}
int main() {
int t = read(), n, m;
init();
while (t--) {
n = read(), m = read();
if (n > m)swap(n, m);
int l = 1, r;
ll ans = 0;
while (l <= n)r = min(n / (n / l), m / (m / l)), ans += (mu[r] - mu[l - 1]) * s[n / l] * s[m / l], l = r + 1;
printf("%lld\n", ans);
}
return 0;
}
Problem b
其实就是第一个问题的修改版，它的下界不再为1。这里需要用到容斥原理。
定义：
定义我们需要求的集合为A，那么有关系：
所以：
后者用一步容斥原理：
所以答案即为：
然后就转化为第一个问题了，直接求四遍即可。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
using namespace std;
typedef long long ll;
int mu[N + 1], prim[N + 1], mark[N + 1], tot, k, a, b, c, d;
inline int read() {
char e = getchar();
int s = 0;
while (e < '-')e = getchar();
while (e > '-')s = (s << 1) + (s << 3) + (e & 15), e = getchar();
return s;
}
inline void getMu() {//求莫比乌斯函数
for (int i = 2; i <= N; i++) {
if (!mark[i])prim[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot; j++) {
if ((ll) i * prim[j] > N)break;
mark[i * prim[j]] = 1;
if (i % prim[j] == 0)break;
else mu[i * prim[j]] = -mu[i];
}
}
mu[1] = 1;
}
ll getAns(int n, int m) {
if (n > m)swap(n, m);
int l = 1, r;
ll ans = 0;
while (l <= n) {
r = min(n / (n / l), m / (m / l));
ans += (ll) (n / (l * k)) * (m / (l * k)) * (mu[r] - mu[l - 1]);
l = r + 1;
}
return ans;
}
int main() {
int t = read();
getMu();
for (int i = 1; i <= N; i++)mu[i] += mu[i - 1];
while (t--) {
a = read(), b = read(), c = read(), d = read(), k = read();
printf("%lld\n", getAns(b, d) - getAns(b, c - 1) - getAns(a - 1, d) + getAns(a - 1, c - 1));
}
return 0;
}
Crash的数字表格
本质上是求：
就是：
枚举最大公因数：
只枚举d的倍数：
利用gcd的第二个套路，得到：
枚举x，这一步转化已经用了很多次了：
只枚举x的倍数，这一步也很常用：
这样不好看，来一步转化：
后面是一个等差数列求和，而$x^2\mu(x)$可以线性筛出，套一个整除分块即可。代码常数大一些，可以卡常过去。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
using namespace std;
typedef long long ll;
int mu[N + 1], prim[N + 1], mark[N + 1], tot, n, m;
inline void getMu() {//求莫比乌斯函数
for (int i = 2; i <= N; i++) {
if (!mark[i])prim[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot; j++) {
if ((ll) i * prim[j] > N)break;
mark[i * prim[j]] = 1;
if (i % prim[j] == 0)break;
else mu[i * prim[j]] = -mu[i];
}
}
mu[1] = 1;
for (int i = 1; i <= N; i++)mu[i] = 1ll * i * i % MOD * 1ll * mu[i] * 1ll % MOD;
for (int i = 1; i <= N; i++)mu[i] = (mu[i] + mu[i - 1]) % MOD;
}
int main() {
cin >> n >> m;
getMu();
if (n > m)swap(n, m);
long long ans = 0ll, inv = 10050505ll;//这是逆元
for (int i = 1; i <= n; i++) {
int p = n / i, l = 1, r;
long long sum = 0, tmp = 1l;
while (l <= p) {
tmp = 1l;
r = min(p / (p / l), (m / i) / (m / i / l));
tmp = tmp * (1ll * mu[r] - 1ll * mu[l - 1]) % MOD;
tmp = tmp * (1ll + p / l) % MOD * 1ll * (p / l) % MOD * 1ll * inv % MOD;
tmp = tmp * (1ll + m / i / l) % MOD * 1ll * (m / i / l) % MOD * 1ll * inv % MOD;
sum = (sum + tmp) % MOD;
l = r + 1;
}
ans = (ans + sum * 1ll * i) % MOD;
}
cout << (ans + MOD) % MOD;
return 0;
}
[NOI2011]能量采集
其实就是求：
重点当然是求前面那一部分，不妨设$m\leq n$。枚举最大公因数：
提出公因子：
然后套路一波：
枚举k（如果你坚持看到了这里，我相信你会想到这一步的）：
其实到这一步已经可以做了，但是这个式子仍然可以优化（本题的亮点），枚举$dk$：
后边这个东西就是$\phi(i)$：
然后就可以$O(n)$地做这题了。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
using namespace std;
int mark[N], prim[N], tot = 0, n, m;
long long euler[N];
int main() {
euler[1] = 1;
for (int i = 2; i <= 100000; i++) {
if (!mark[i])prim[++tot] = i, euler[i] = i - 1;
for (int j = 1; j <= tot; j++) {
if (i * prim[j] > 100000)break;
mark[i * prim[j]] = 1;
if (i % prim[j] == 0) {
euler[i * prim[j]] = euler[i] * prim[j];
break;
} else euler[i * prim[j]] = euler[i] * (prim[j] - 1);
}
}
for (int i = 2; i <= 100000; i++)euler[i] += euler[i - 1];
cin >> n >> m;
if (m > n)swap(n, m);
int l = 1, r;
long long ans = 0;
while (l <= m) {
r = min(m / (m / l), n / (n / l)), ans += 1ll * (n / l) * (m / l) * (euler[r] - euler[l - 1]), l = r + 1;
}
cout << 2ll * ans - 1ll * n * m;
return 0;
}
简单的数学题
简单你妹夫，就是求：
枚举一下最大公因数：
右边的式子拿出来，套路一波。设：
那么：
所以：
其实$f(d)$是可以直接求的，它就是：
记$S(x)=1+2+\cdots+x$，那么就是：
代入原式，得到：
考虑每一个倍数，因子对它的贡献，可以得到：
后面那个东西就是莫比乌斯函数与单位函数的狄利克雷卷积，就是欧拉函数：
前面那一个当然可以整除分块，现在问题就是求后面的前缀和，但是$n$高达$10^{10}$，用线性筛会T，于是考虑杜教筛。
写出杜教筛的套路式：
令$g(x)=x^2$，那么给$g(x)$和$x^2\phi(x)$做狄利克雷卷积：
然后：
这些$i^3$和$d^2$的求和式都是能$O(1)$求的，这样就可以用杜教筛求前缀和，再套上上面的整除分块，本题就解决了。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
using namespace std;
unordered_map<long long, long long> mmp;
long long n, prim[N + 5], tot = 0, ps[N + 5], p, inv6, inv4;
bool vis[N + 5];
long long inv(int x) {
long long ans = 1, sta = x, y = p - 2;
while (y) {
if (y & 1)ans = ans * sta % p;
sta = sta * sta % p, y >>= 1;
}
return ans;
}
inline void init() {
ps[1] = 1;
for (int i = 2; i <= N; i++) {
if (!vis[i])ps[i] = i - 1, prim[++tot] = i;
for (int j = 1; j <= tot; j++) {
if (prim[j] * i > N)break;
vis[i * prim[j]] = true;
if (i % prim[j]) ps[i * prim[j]] = ps[i] * (prim[j] - 1) % p;
else {
ps[i * prim[j]] = ps[i] * prim[j] % p;
break;
}
}
}
for (int i = 1; i <= N; i++)ps[i] = 1ll * i * i % p * ps[i] % p;
for (int i = 2; i <= N; i++)ps[i] = (ps[i] + ps[i - 1]) % p;
}
inline long long S(long long x) {
x %= p;
return x * (x + 1) % p * (2 * x + 1) % p * inv6 % p;
}
inline long long T(long long x) {
x %= p;
return x * x % p * (x + 1) % p * (x + 1) % p * inv4 % p;
}
long long get(long long x) {
if (x <= N)return ps[x];
if (mmp.count(x) != 0)return mmp[x];
long long ans = 0, l = 2, r;
while (l <= x)r = x / (x / l), ans += (S(r) - S(l - 1)) * get(x / l) % p, ans %= p, l = r + 1;
return mmp[x] = (T(x) - ans) % p;
}
int main() {
cin >> p >> n;
init(), inv6 = inv(6), inv4 = inv(4);
long long ans = 0, l = 1, r;
while (l <= n) {
r = n / (n / l);
ans = (ans + T(n / l) * (get(r) - get(l - 1)) % p) % p, l = r + 1;
}
cout << (ans + p) % p;
return 0;
}
[CQOI2015]选数
设在其中选$n$个数$gcd=d$的方案数为$f(d)$，$gcd=kd$的方案数为$g(d)$，那么有：
然后，喜闻乐见的一步：
$g(d)$能直接求：
代入原式：
枚举倍数：
后面要做什么就很显然了，这里的前缀和需要快速求，得用杜教筛。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
using namespace std;
map<int, int> mmp;
int mu[N + 1], mark[N + 1], prim[N + 1], tot = 0;
int n, k, L, H;
inline int qPow(long long x, int y) {
long long ans = 1, sta = x;
while (y) {
if (y & 1)ans = ans * sta % MOD;
sta = sta * sta % MOD, y >>= 1;
}
return ans;
}
inline void init() {
mu[1] = 1;
for (int i = 2; i <= N; i++) {
if (!mark[i])prim[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot; j++) {
if (1ll * i * prim[j] > N)break;
mark[i * prim[j]] = 1;
if (i % prim[j] == 0)break;
else mu[i * prim[j]] = -mu[i];
}
}
for (int i = 2; i <= N; i++)mu[i] += mu[i - 1];
}
int get(long long x) {//杜教筛：莫比乌斯函数
if (x <= N)return mu[x];
if (mmp.count(x) != 0)return mmp[x];
long long ans = 0, l = 2, r;
while (l <= x)r = x / (x / l), ans += (r - l + 1) * get(x / l), l = r + 1;
return mmp[x] = 1 - ans;
}
int main() {
cin >> n >> k >> L >> H;
init();
int sj = H / k, l = 1, r;
long long ans = 0;
while (l <= sj) {
if ((L - 1) / l != 0)r = min(H / (H / l), (L - 1) / ((L - 1) / l));
else r = H / (H / l);
ans = (ans + (get(r) - get(l - 1)) * qPow(H / (1ll * l * k) - (L - 1) / (1ll * l * k), n) % MOD) % MOD;
l = r + 1;
}
cout << (ans + MOD) % MOD;
return 0;
}