0%

浅谈对莫比乌斯反演的理解

莫比乌斯反演通常用于解决求和内含$\gcd$的问题,通常是多个求和公式嵌套而成并且需要在线性的复杂度(甚至需要低于线性时间复杂度)内解决该问题。该博文就是来谈谈对莫比乌斯反演的理解。

基本内容

莫比乌斯函数

莫比乌斯函数是指以下的函数:
$$\mu(n)=\begin{cases}1,&n=1\\(-1)^k,&\text{若}n\text{无平方数因数,且}n=p_1p_2p_3…p_n\\0,&\text{若}n\text{有大于}1\text{的平方数因数}\end{cases}$$

同时,莫比乌斯函数是一个积性函数

$$\sum\limits_{d|n}\mu(d)=\begin{cases}1,&\text{若}n=1\\0,&\text{其他情况}\end{cases}=[n=1]$$

莫比乌斯反演公式

假设对于数论函数$f(n)$和$F(n)$,有以下关系式:
$$F(n)=\sum\limits_{d|n}f(d)$$
则其莫比乌斯反演公式定义为:
$$f(n)=\sum\limits_{d|n}\mu(d)F(\frac{n}{d})$$
它还有第二种形式:
$$F(n)=\sum\limits_{n\mid d}f(d)\Rightarrow f(n)=\sum\limits_{n\mid d}\mu(\frac{d}{n})F(d)$$

杜教筛

杜教筛是用来简化求积性函数前缀和的一种方法,也用到了下面所说的黑科技。杜教筛不是本文的重点,所以就不多介绍了。(学习杜教筛可以到这篇博文,非常详细)主要是我不会讲。

举个例子

莫比乌斯反演可以解决一类求和中含有求最大公因数(以下简称$\gcd$)的问题。
举一个稍简单的例子,让你求这样的一个式子:
$$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gcd(i,j)\text{ ①}$$
$1\le n\le10^{9}$。

显然暴力不可行,那我们试试推导一番。

推导公式

首先枚举所有$\gcd$得到的值:
$$\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[\gcd(i,j)=d]\text{ ②}$$

我们知道,$\gcd(i,j)=d$时,$i$和$j$一定是$d$的倍数,那么我们在已知$d$的前提下,$i$和$j$都只需要枚举$d$的倍数即可,即上式可化为:

$$\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}[\gcd(i,j)=1]\text{ ③}$$

显然右边的式子是比较整齐的,令$d$右边的求和式整理为一个函数:

$$f(m,k)=\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{m}[\gcd(i,j)=k]\text{ ④}$$

为什么要把③中的$\gcd(i,j)=1$中的$1$也变成一个变量处理呢?这是因为处理这个函数需要使用莫比乌斯反演公式,而把这个$1$作为变量更好理解。

既然我们需要用莫比乌斯反演公式,那是否能找到一个函数使得:

$$F(m,k)=\sum\limits_{k|l}f(m,l)\text{ ⑤}$$

您好,有的。

$$F(m,k)=\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{m}[\gcd(i,j)=k\text{的倍数}]=(\lfloor\frac{m}{k}\rfloor)^2\text{ ⑥}$$

这一步非常难以理解,请读者耐心思考。

第一个难以理解的点在于⑤式中如何求到所有的$l$能整数$k$。仔细分析就知道这个$l$有无穷个,但在无穷个合适的取值中大部分的$l$代入函数f中的值都为0。
第二个难以理解的点在于⑥式中为什么函数F是这样的真的会使得⑤式成立吗。因为所有$l$组成的集合就是所有$k$的倍数的集合。

这个函数F的值是非常容易求的,那么得到函数F之后,那么就可以利用莫比乌斯反演公式迎刃而解了。

$$\text{原式}=\sum\limits_{d=1}^{n}d\cdot f(\lfloor\frac{n}{d}\rfloor,1)=\sum\limits_{d=1}^{n}d\sum\limits_{1|i}\mu(i)F(\lfloor\frac{n}{d}\rfloor,i)$$

$$=\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{\infty}\mu(i)F(\lfloor\frac{n}{d}\rfloor,i)=\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{\infty}\mu(i)(\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{i}\rfloor)^2$$

$$=\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\mu(i)(\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{i}\rfloor)^2\text{ ⑦}$$

黑科技

显然⑦式中$\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}…$的部分只与$\lfloor\frac{n}{d}\rfloor$的值有关。
当$d\le\sqrt{n}$时,$\lfloor\frac{n}{d}\rfloor$显然只有$O(\sqrt{n})$个取值;当$d\gt n$时,$\lfloor\frac{n}{d}\rfloor\lt\sqrt{n}$显然也只有$O(\sqrt{n})$个取值;对于固定的$\lfloor\frac{n}{d}\rfloor$,$d$的取值是一段连续的区间,这段区间是$[\left\lfloor\frac{n}{\left\lfloor\frac{n}{d}\right\rfloor+1}\right\rfloor+1,\left\lfloor\frac{n}{\left\lfloor\frac{n}{d}\right\rfloor}\right\rfloor]$,因此可以$O(\sqrt{n})$计算所求。
即设$f(n,i)=\lfloor\frac{n}{i}\rfloor$,那么使得$f(n,i)=f(n,j)$的$j\in[\left\lfloor\frac{n}{\left\lfloor\frac{n}{i}\right\rfloor+1}\right\rfloor+1,\left\lfloor\frac{n}{\left\lfloor\frac{n}{i}\right\rfloor}\right\rfloor]$。
还有一种是向上取整的,即设$f(n,i)=\lceil\frac{n}{i}\rceil$,那么使得$f(n,i)=f(n,j)$的$j\in[\lceil\frac{n}{\lceil\frac{n}{i}\rceil}\rceil,\left\lceil\frac{n}{\left\lceil\frac{n}{i}\right\rceil-1}\right\rceil-1]$。(若$\lceil\frac{n}{i}\rceil=1$,上界为无穷大)

那么上面的式子还有一个难点,就是求莫比乌斯函数的前缀和,用杜教筛就可以求出来。最终总的复杂度是$O(n^{\frac{3}{4}})$。如果$1\le n\le 10^{10}$,需要使用其他方法化简上面的式子变成欧拉函数前缀和同样也用杜教筛求。

做个小结

通常情况下:
第一步都需要枚举所有可能的$\gcd$的值
第二步化简成具有向上取整和向下取整的公式;
第三步利用莫比乌斯反演来使公式更容易求出;
第四步看是否能以及需要使用杜教筛来简化运算。

可以说这种题目非常地套路了。

例题

ECNU3308 lcm 与 gcd (1)

题意:给你下面的式子,让你去求它的值:

$$\sum_{i=1}^n \sum_{j=1}^m \frac{\mathrm{lcm}(i,j)}{\gcd(i,j)} \bmod 1000000~007$$

尝试着去推一下吧。详细解析

牛牛的最大公约数

题意:给你下面的式子,让你去求它的值:

$$\sum\limits_{i_1=L}^{R}\sum\limits_{i_2=L}^{R}\sum\limits_{i_3=L}^{R}…\sum\limits_{i_N=L}^{R}[gcd(i_1,i_2,i_3,…,i_n)=K]$$

化简:

第二步:

$$f(d)=\sum\limits_{i_1=\lceil\frac{L}{K}\rceil}^{\lfloor\frac{R}{K}\rfloor}\sum\limits_{i_2=\lceil\frac{L}{K}\rceil}^{\lfloor\frac{R}{K}\rfloor}…\sum\limits_{i_N=\lceil\frac{L}{K}\rceil}^{\lfloor\frac{R}{K}\rfloor}[gcd(i_1,i_2,…,i_N)=d]$$

即答案为f(1)。

第三步:

$$F(d)=\sum\limits_{i_1=\lceil\frac{L}{K}\rceil}^{\lfloor\frac{R}{K}\rfloor}\sum\limits_{i_2=\lceil\frac{L}{K}\rceil}^{\lfloor\frac{R}{K}\rfloor}…\sum\limits_{i_N=\lceil\frac{L}{K}\rceil}^{\lfloor\frac{R}{K}\rfloor}[gcd(i_1,i_2,…,i_N)=d的倍数]=(\lfloor\frac{\lfloor\frac{R}{K}\rfloor-\lceil\frac{L}{K}\rceil+1}{d}\rfloor)^N$$

则:

$f(1)=\sum\limits_{d=1}^{\infty}\mu(i)\cdot(\lfloor\frac{\lfloor\frac{R}{K}\rfloor-\lceil\frac{L}{K}\rceil+1}{d}\rfloor)^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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#include<bits/stdc++.h>
using namespace std;
long long fpow(long long a, long long b, const long long &p) {
long long ans = 1ll;
for (a %= p; b; a = a * a % p, b >>= 1) {
if (b & 1) ans = ans * a % p;
}
return ans;
}
const long long MOD = 1e9 + 7;
const int MAXN = 1e7 + 5;
bool isnpri[MAXN];
vector<int>pri;
int mu[MAXN];
void mobius() {
mu[1] = 1;
for(int i = 2; i < MAXN; i++) {
if(!isnpri[i]) {
pri.push_back(i);
mu[i] = -1;
}
for(int j = 0; j < signed(pri.size()); j++) {
const int cur = i * pri[j];
if(cur >= MAXN) {
break;
}
isnpri[cur] = true;
if(i % pri[j] == 0) {
mu[cur] = 0;
break;
} else {
mu[cur] = -mu[i];
}
}
}
for(int i = 2; i < MAXN; i++) {
mu[i] += mu[i - 1];
}
}
unordered_map<int, int> mp;
int pm(int x) {
if(x < MAXN) {
return mu[x];
}
unordered_map<int, int>::iterator it = mp.find(x);
if(it != mp.end()) {
return it->second;
}
long long ans = 1;
for(int i = 2, j; i <= x; i = j + 1) {
j = x / (x / i);
ans = (((ans - 1ll * pm(x / i) * (j - i + 1)) % MOD) + MOD) % MOD;
}
mp[x] = ans;
return ans;
}
int main() {
mobius();
int n, k, l, r;
cin >> n >> k >> l >> r;
long long ans = 0;
const int A = r / k, B = (l + k - 1) / k;
for(int d = 1, p; d <= A; d = p + 1) {
const int C = (B + d - 1) / d;
const int D = C == 1 ? INT_MAX : (B + C - 2) / (C - 1) - 1;
p = min(A / (A / d), D); // 求最大的d使某部分相同,而另一部分通过求函数前缀和获得。
const long long cur = A / d - C + 1;
if(cur <= 0) continue;
ans = (((ans + fpow(cur, n, MOD) * (pm(p) - pm(d - 1))) % MOD) + MOD) % MOD;
}
cout << ans << endl;
}

参考资料

作者:ConanYu
原文地址:https://conanyu.github.io/2019/04/29/mobius/

回到开头