莫比乌斯反演

前置知识

莫比乌斯函数

定义

设正整数

n = p_1^{a_1} p_2^{a_2} \dots p_k^{a_k}

其中

p_i

是互不相同的质数，

a_i \ge 1

。

则莫比乌斯函数

\mu(n)

定义为：

\mu(n) = \begin{cases} 1 & \text{若 } n = 1 \\ (-1)^k & \text{若 } a_1 = a_2 = \dots = a_k = 1 \\ 0 & \text{若 } \exists \, a_i > 1 \end{cases}

性质

这里我们介绍最重要的两个性质：

莫比乌斯函数为积性函数，即如果 $a \perp b$ ，则 $\mu(ab)=\mu(a)\mu(b)$

这里我们使用最基础的方法来证明一下。

假设

\gcd(a, b) = 1

。

若 $a$ 或 $b$ 中有一个包含平方因子。则 $ab$ 必然也包含平方因子。此时 $\mu(a)\mu(b) = 0$ ，而 $\mu(ab) = 0$ 。等式成立。
若 $a$ 和 $b$ 都不含平方因子。设 $a = p_1 p_2 \dots p_k$ ， $b = q_1 q_2 \dots q_m$ 。因为 $\gcd(a, b) = 1$ ，所以所有的质数 $\{p_i\}$ 和 $\{q_j\}$ 互不相同。那么 $ab$ 就是 $k+m$ 个不同质数的乘积。根据定义：

\mu(a) = (-1)^k, \quad \mu(b) = (-1)^m

\mu(ab) = (-1)^{k+m} = (-1)^k \cdot (-1)^m = \mu(a)\mu(b)

证明完毕。

还有一种用狄利克雷卷积证明的办法，感兴趣的读者可以下去自行探索。

这个性质解释了为什么我们能够用线性筛把

1\sim n

之间的

\mu

值在

O(n)

的复杂度内算出。

具体的，当我们处理

i \times prime[j]

时，

如果两者互质，利用积性性质： $mu[i \times prime[j]] = mu[i] \times mu[prime[j]]$ ，

因为

prime[j]

是质数，

\mu(prime[j])

永远等于

-1

。

所以

\mu(i \times prime[j]) = \mu(i) \times (-1) = -\mu(i)

。

如果 $prime[j] \mid i$ ，说明出现了平方因子 $prime[j]^2$ ，根据定义，结果直接为 $0$ 。

线性筛代码

CPP

inline void init(int limit)
{
    fill(isprime + 1, isprime + 1 + limit, true);

    miu[1] = 1;

    for (int i = 2; i <= limit; ++i)
    {
        if (isprime[i]) { prime[++sub] = i; miu[i] = -1; }

        for (int j = 1; j <= sub && i * prime[j] <= limit; ++j)
        {
            isprime[i * prime[j]] = false;

            if (i % prime[j] == 0)
            {
                miu[i * prime[j]] = 0;
                break;
            }
            else miu[i * prime[j]] = -miu[i];
        }
    }
}

因数和性质

\sum_{d|n} \mu(d) = [n=1] = \begin{cases} 1 & \text{若 } n = 1 \\ 0 & \text{若 } n > 1 \end{cases}

这个性质在莫比乌斯反演中非常的重要。

比如，当我们遇到“当且仅当

\gcd(i, j) = 1

时贡献 1，否则贡献 0”的情况。我们发现判定

\gcd(i, j)=1

非常的困难，这时候我们可以转换为

[\gcd(i, j)=1] = \sum_{d|\gcd(i, j)} \mu(d)

然后我们就可以通过交换求和次序，把复杂的

\gcd

计数转化为简单的倍数统计，这是所有莫比乌斯反演题目的通用套路。

所以，当我们遇到

\gcd(i, j)=d

的问题，这道题有很大的概率要运用莫比乌斯反演进行解决。

证明的过程也比较简洁。设

n = p_1^{a_1} \dots p_k^{a_k}

，其中

k

为不同质因子的个数。

由于当因数

d

含有平方因子时

\mu(d)=0

，因此有效的

d

只能是这

k

个质因子的某种组合。

我们将每一项

\mu(d)

的贡献看作是从

k

个质因子中选取

r

个的结果：

\sum_{d|n} \mu(d) = \sum_{r=0}^{k} \binom{k}{r} (-1)^r

根据二项式定理

(x+y)^k = \sum_{r=0}^k \binom{k}{r} x^{k-r} y^r

，我们令

x=1, y=-1

：

\sum_{r=0}^{k} \binom{k}{r} (-1)^r = (1 - 1)^k = \begin{cases} 1 & k=0 \text{ (即 } n=1) \\ 0 & k \ge 1 \text{ (即 } n>1) \end{cases}

还有一些其他有意思的性质，由于做题的时候用不到，这里我们就不再介绍了。

定义

我们做题的时候一般把我们要求的目标函数定义为

f(d)

，再定义一个辅助函数

F(n)

，一般这个函数都比较好计算，然后

F(n)

和

f(d)

之间存在某种关系，最后通过

F(n)

反推出

f(d)

，然后通过一系列操作降低复杂度最后求得结果。

1. 约数反演

若有

F(n) = \sum_{d|n} f(d)

，则有：

f(n) = \sum_{d|n} \mu(d) F\left(\frac{n}{d}\right)

证明：

首先给大家介绍几个东西。

狄利克雷卷积定义： $(a * b)(n) = \sum_{d|n} a(d)b(\frac{n}{d})$ 。它有着代数的运算性质
单位元： $\epsilon(n) = [n=1]$ 。任何函数 $f * \epsilon = f$ 。
恒等函数： $1(n) = 1$ 。

已知

F = f * 1

（

F(n) = \sum_{d|n} f(d)

）。

F * \mu = (f * 1) * \mu

F * \mu = f * (1 * \mu)

容易得知，

\mu * 1 = \epsilon

。

F * \mu = f * \epsilon = f

f = F * \mu

，即

f(n) = \sum_{d|n} \mu(d) F(\frac{n}{d})

。证毕。

2. 倍数反演

倍数反演是竞赛中最常用的反演，绝大多数问题中运用的都是此反演。

若有

F(n) = \sum_{n|d} f(d)

，则有：

f(n) = \sum_{n|d} \mu\left(\frac{d}{n}\right) F(d)

证明过程与上面相似，这里不再阐述。

3.指数反演

这个基本上不会用，貌似数学竞赛中有时候会用到。

若有

F(n) = \prod_{d|n} f(d)

，则有：

f(n) = \prod_{d|n} F\left(\frac{n}{d}\right)^{\mu(d)}

做题思路

一般的思路是先将一个问题抽象成

\gcd

计数模型，

得到一个类似这样的式子：

Ans = \sum_{i=1}^{n} \sum_{j=1}^{m} [\gcd(i, j) = k]

令

i' = i/k, j' = j/k

Ans = \sum_{i'=1}^{\lfloor n/k \rfloor} \sum_{j'=1}^{\lfloor m/k \rfloor} [\gcd(i', j') = 1]

令

N = \lfloor n/k \rfloor, M = \lfloor m/k \rfloor

。

设

f(d)

：表示

\gcd(i, j)

恰好等于

d

的对数（我们最终想要

f(1)

）。

F(d)

：表示

\gcd(i, j)

是 $d$ 的倍数 的对数。

显然有

F(d) = \sum_{d|k} f(k)

又可以推出

F(d) = \lfloor \frac{N}{d} \rfloor \cdot \lfloor \frac{M}{d} \rfloor

f(1) = \sum_{1|d} \mu(\frac{d}{1}) F(d) \implies \sum_{d=1}^{\min(N,M)} \mu(d) \cdot \lfloor \frac{N}{d} \rfloor \cdot \lfloor \frac{M}{d} \rfloor

然后我们通过预处理或者整除分块操作来降低复杂度完成。

有时候我们需要用到求和顺序交换，提出因数

d

，结合狄利克雷卷积等经典操作，具体的在后面题目中会讲解。

简单的变形，如求

\sum \sum \gcd(i, j)

，我们可以枚举

\gcd

的值：

\sum_{g=1}^{\min(n, m)} g \cdot \sum_{i=1}^n \sum_{j=1}^m [\gcd(i, j) = g]

。后面就一直推下去就好。

整除分块

如果我们要计算

Ans = \sum_{i=1}^{n} \lfloor \frac{n}{i} \rfloor

我们可以发现

\lfloor n/i \rfloor

的值在一段连续的区间内是不变的。

我们想到可以把这些相同的块一次性算出来，就可以降低一定的复杂度。

假设当前块的左端点是

l

，我们需要找到最大的右端点

r

，使得对于所有

i \in [l, r]

，

\lfloor n/i \rfloor

都相等。这个

r

的计算公式极其简洁：

r = \lfloor \frac{n}{\lfloor n/l \rfloor} \rfloor

证明（简述）：设

k = \lfloor n/l \rfloor

。我们要找最大的

r

满足

r \cdot k \le n

，显然

r \le n/k

。向下取整即得

r = \lfloor n/k \rfloor

。

参考代码

CPP


    long long ans = 0;

    for (int l = 1, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        // 当前块的值是 n/l，块的长度是 (r - l + 1)
        ans += (long long)(r - l + 1) * (n / l);
    }

下面讲解一些莫比乌斯反演好题，让大家能对这一算法有着更深的理解。

$\textit{\textbf {Eg1}}$

题目描述：P2522 [HAOI2011] Problem b

这道题是 莫比乌斯反演 的经典入门题，所以拿出来作为讲解的第一道题。

题目简述

题目要求计算满足

a \le x \le b

，

c \le y \le d

且

\gcd(x, y) = k

的数对

(x, y)

的个数。

推导过程

由于区间

[a, b]

和

[c, d]

不是从 1 开始的，直接计算比较麻烦。我们可以利用二维前缀和的思想（容斥原理）将其转化为四个从

(1, 1)

开始的子问题。

定义函数

solve(n, m, k)

为：满足

1 \le x \le n

，

1 \le y \le m

且

\gcd(x, y) = k

的数对数量。那么原问题的答案为：

Ans = \text{solve}(b, d, k) - \text{solve}(a-1, d, k) - \text{solve}(b, c-1, k) + \text{solve}(a-1, c-1, k)

现在我们的核心任务是快速求出

solve(n, m, k)

。

\text{solve}(n, m, k) = \sum_{x=1}^{n} \sum_{y=1}^{m} [\gcd(x, y) = k]

\gcd(x, y) = k

等价于

\gcd(x/k, y/k) = 1

。令

n' = \lfloor n/k \rfloor, m' = \lfloor m/k \rfloor

。问题转化为求

1 \le i \le n', 1 \le j \le m'

且

\gcd(i, j) = 1

的对数。

\text{solve}(n, m, k) = \sum_{i=1}^{n'} \sum_{j=1}^{m'} [\gcd(i, j) = 1]

我们定义

$f(d)$ ：满足 $\gcd(i, j) = d$ 的对数。
$F(d)$ ：满足 $d \mid \gcd(i, j)$ 的对数

显然有关系：

F(d) = \sum_{d|K} f(K)

。

同时，

d \mid \gcd(i, j)

意味着

i

是

d

的倍数且

j

是

d

的倍数。在

1 \dots n'

中，

d

的倍数有

\lfloor n'/d \rfloor

个；在

1 \dots m'

中有

\lfloor m'/d \rfloor

个。故直接得出：

F(d) = \lfloor \frac{n'}{d} \rfloor \lfloor \frac{m'}{d} \rfloor

根据倍数形式的反演公式：

f(d) = \sum_{d|K} \mu(\frac{K}{d}) F(K)

。取

d=1

，代入

F(K)

：

f(1) = \sum_{K=1}^{\min(n', m')} \mu(K) \cdot \lfloor \frac{n'}{K} \rfloor \lfloor \frac{m'}{K} \rfloor

虽然推导出了

f(1) = \sum_{i=1}^{\min(n', m')} \mu(i) \lfloor \frac{n'}{i} \rfloor \lfloor \frac{m'}{i} \rfloor

，但单次询问

O(N)

会超时，我们考虑用怎么优化。

我们可以发现后面那一堆是整除分块的经典形式，即对于

\lfloor n'/i \rfloor

，当

i

在一定范围内变化时，其值是不变的。

我们可以预处理

\mu

的前缀和

S(i) = \sum_{j=1}^i \mu(j)

。对于每一块

[l, r]

，其贡献为：

(S(r) - S(l-1)) \cdot \lfloor \frac{n'}{l} \rfloor \cdot \lfloor \frac{m'}{l} \rfloor

其中右边界

r = \min(n' / (n'/l), m' / (m'/l))

。

最终

solve(n, m, k)

的计算公式为：

solve(n, m, k) = \sum_{i=1}^{\min(\lfloor n/k \rfloor, \lfloor m/k \rfloor)} \mu(i) \lfloor \frac{n}{ki} \rfloor \lfloor \frac{m}{ki} \rfloor

参考代码

CPP

#include <bits/stdc++.h>

#define ll long long

using namespace std;

const int MAX = 5e4 + 10;

int Q;

int prime[MAX], sub;
bool isprime[MAX];
int miu[MAX], sum[MAX];

inline void init(int limit)
{
    fill(isprime + 1, isprime + 1 + limit, true);

    isprime[1] = false;
    miu[1] = 1;

    for (int i = 2; i <= limit; ++i)
    {
        if (isprime[i]) { prime[++sub] = i; miu[i] = -1; }

        for (int j = 1; j <= sub && i * prime[j] <= limit; ++j)
        {
            isprime[i * prime[j]] = false;

            if (i % prime[j] == 0)
            {
                miu[i * prime[j]] = 0;
                break;
            }
            else miu[i * prime[j]] = -miu[i];
        }
    }

    for (int i = 1; i <= limit; ++i) sum[i] = sum[i - 1] + miu[i];
}

inline ll calc(int n, int m)
{
    if (n > m) swap(n, m);
    
    ll res = 0;

    for (int l = 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        res += (ll)(sum[r] - sum[l - 1]) * (n / l) * (m / l);
    }

    return res;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0); cout.tie(0);

    init(50000);

    cin >> Q;

    while (Q--)
    {
        int a, b, c, d, k; cin >> a >> b >> c >> d >> k;

        // 核心逻辑：容斥原理，利用 calc 计算
        // 记得除以 k 转化成 gcd(x,y)=1 的问题
        ll ans = calc(b / k, d / k) 
               - calc((a - 1) / k, d / k) 
               - calc(b / k, (c - 1) / k) 
               + calc((a - 1) / k, (c - 1) / k);

        cout << ans << "\n";
    }

    return 0;
}

$\textit{\textbf {Eg2}}$

题目描述：P2158 [SDOI2008] 仪仗队

这道题是一道经典的数论题目。要解决它，我们需要将视觉问题转化为数学模型。

推导过程

1. 建立坐标系

假设 C 君的位置在坐标原点 $(0,0)$ 。那么学生方阵占据的坐标范围是从 $(0,0)$ 到 $(N-1, N-1)$ 。

2. 判定“能看到”的条件

一个点

(x, y)

能被看到，当且仅当在 C 君的视线上，它是第一个点。如果存在另一个点

(x', y')

在连线

O(x, y)

上且更靠近原点，那么

(x, y)

就会被挡住。这意味着

x

和

y

必须互质。

如果 $\gcd(x, y) = d > 1$ ，那么点 $(\frac{x}{d}, \frac{y}{d})$ 就会挡住 $(x, y)$ 。
因此，问题转化为：在 $0 \le x, y < N$ 的范围内，有多少对 $(x, y)$ 满足 $\gcd(x, y) = 1$ ？

3. 特殊情况处理

当 $N=1$ 时，方阵里其实没人，或者说只有 C 君自己，输出 $0$ （或者根据具体定义可能是 $0$ ）。但在本题数据范围内， $N=1$ 需要特判，通常这种视野题 $N=1$ 结果为 $0$ 。
点 $(1,0)$ ， $(0,1)$ 和 $(1,1)$ 是必然能看到的。

推导式子

这道题能单纯用欧拉函数推导出来，但是为了练习莫比乌斯反演，这里我们从莫比乌斯反演的角度来推导一下。

我们先不考虑坐标轴上的点。

设

n = N-1

，我们需要求：

\sum_{i=1}^{n} \sum_{j=1}^{n} [\gcd(i, j) = 1]

定义 $f(d)$ ：满足 $\gcd(x, y) = d$ 的数对 $(x, y)$ 的个数，其中 $1 \le x, y \le n$ 。
定义 $F(d)$ ：满足 $d | \gcd(x, y)$ 的数对 $(x, y)$ 的个数，其中 $1 \le x, y \le n$ 。

显然，如果一个数对的公约数是

d

的倍数，那么它们的最大公约数一定是

d, 2d, 3d \dots

中的某一个。

F(d) = \sum_{d|k, k \le n} f(k)

同时，根据定义，

1 \dots n

范围内

d

的倍数有

\lfloor \frac{n}{d} \rfloor

个。要让

x

和

y

都是

d

的倍数，组合数就是：

F(d) = \lfloor \frac{n}{d} \rfloor^2

根据莫比乌斯反演的第二种形式（倍数形式）：若

F(d) = \sum_{d|k} f(k)

，则有：

f(d) = \sum_{d|k, k \le n} \mu(\frac{k}{d}) F(k)

我们要求的是

f(1)

，代入公式：

f(1) = \sum_{1|k, k \le n} \mu(\frac{k}{1}) F(k)

将

F(k) = \lfloor \frac{n}{k} \rfloor^2

代入上式，得到最终计算公式：

f(1) = \sum_{k=1}^{n} \mu(k) \lfloor \frac{n}{k} \rfloor^2

说明

最后只需要加上坐标轴上的两个点即可

参考代码

CPP

#include <bits/stdc++.h>

using namespace std;

const int N = 4e4 + 10;
const int MAX = 4e4;

int n;

bool isprime[N];
int miu[N], prime[N], sub;

inline void initmiu()
{
    memset(isprime, true, sizeof isprime);

    miu[1] = 1;

    for (int i = 2; i <= MAX; ++i)
    {
        if (isprime[i]) { miu[i] = -1; prime[++sub] = i; }

        for (int j = 1; j <= sub && i * prime[j] <= MAX; ++j)
        {
            isprime[i * prime[j]] = false;

            if (i % prime[j] == 0) { miu[i * prime[j]] = 0; break; }
            else miu[i * prime[j]] = -miu[i];
        }
    }
}

int main()
{
    scanf("%d", &n);

    if (n == 1) { puts("0"); exit(0); }

    initmiu();

    long long ans = 0;

    n--;

    for (int i = 1; i <= n; ++i)
    {
        ans += (long long) miu[i] * (n / i) * (n / i);
    }

    printf("%lld\n", ans + 2);

    return 0;
}

$\textit{\textbf {Eg3}}$

题目描述：P1829 [国家集训队] Crash的数字表格 / JZPTAB

题目简述

求解

\text{Ans} = \sum_{i=1}^n \sum_{j=1}^m \text{lcm}(i, j) \pmod{20101009}

推导过程

利用基本性质

\text{lcm}(i, j) = \frac{i \cdot j}{\gcd(i, j)}

，公式变为：

\text{Ans} = \sum_{i=1}^n \sum_{j=1}^m \frac{i \cdot j}{\gcd(i, j)}

依照一贯的做法，我们看到

\gcd(i, j)

可以选择提取出来

\gcd(i, j)

的所有取值。

这里我们枚举最大公约数

d

。当

\gcd(i, j) = d

时，这一项的贡献是

\frac{i \cdot j}{d}

。

公式变形为:

\text{Ans} = \sum_{d=1}^{\min(n, m)} \frac{1}{d} \sum_{i=1}^n \sum_{j=1}^m [\gcd(i, j) = d] \cdot i \cdot j

我们再把

d

提取出来，式子变为:

\text{Ans} = \sum_{d=1}^{\min(n, m)}d\sum_{i=1}^{n/d} \sum_{j=1}^{m/d} [\gcd(i, j) = 1] \cdot i \cdot j

看得出来，后面这一堆是莫比乌斯反演的经典应用，因此我们考虑怎么优化

\sum_{i=1}^{n/d} \sum_{j=1}^{m/d} [\gcd(i, j) = 1] \cdot i \cdot j

根据莫比乌斯反演的应用，我们一般设

f(d)

为目标函数，所以这里我们设

f(d)

为

\gcd(i, j)

恰好等于

d

的

i \cdot j

之和：

f(d) = \sum_{i=1}^n \sum_{j=1}^m [\gcd(i, j) = d] \cdot i \cdot j

而我们的目标就是

f(1)

。

按照反演的套路，我们需要定义辅助函数

G(k)

：

G(k) = \sum_{k|d} f(d) = \sum_{i=1}^n \sum_{j=1}^m [k \mid \gcd(i, j)] \cdot i \cdot j

然后我们能发现这里

G(k)

的意义是:

G(k)

为

\gcd(i, j)

是

k

的倍数的

i \cdot j

之和。

那么按照反演的倍数公式，

f(d) = \sum_{d|k} \mu(\frac{k}{d}) G(k)

我们现在考虑如何计算

F(k)

进而计算出我们的目标函数。

条件

k \mid \gcd(i, j)

等价于

k|i

且

k|j

。

我们可以设

i = k \cdot a, j = k \cdot b

，代入式子：

G(k) = \sum_{a=1}^{\lfloor n/k \rfloor} \sum_{b=1}^{\lfloor m/k \rfloor} (k \cdot a) \cdot (k \cdot b)

提取出来

k

：

G(k) = k^2 \left( \sum_{a=1}^{\lfloor n/k \rfloor} a \right) \left( \sum_{b=1}^{\lfloor m/k \rfloor} b \right)

如果我们令等差数列求和函数

g(x) = \frac{x(x+1)}{2}

，那么：

G(k) = k^2 \cdot g(\lfloor n/k \rfloor) \cdot g(\lfloor m/k \rfloor)

此时我们将将

F(k)

代入

f(d)

f(d) = \sum_{d|k} \mu(\frac{k}{d}) \cdot k^2 \cdot g(\lfloor n/k \rfloor) \cdot g(\lfloor m/k \rfloor)

为了方便计算，令

k = d \cdot t

：

f(d) = \sum_{t=1}^{\min(n/d, m/d)} \mu(t) \cdot (dt)^2 \cdot g(\lfloor \frac{n}{dt} \rfloor) \cdot g(\lfloor \frac{m}{dt} \rfloor)

我们最终的答案：

\text{Ans} = \sum_{d=1}^{\min(n, m)} \frac{1}{d} \cdot f(d)

把刚才推导的

f(d)

代入：

\text{Ans} = \sum_{d=1} \frac{1}{d} \sum_{t=1} \mu(t) \cdot d^2 t^2 \cdot g(\lfloor \frac{n}{dt} \rfloor) \cdot g(\lfloor \frac{m}{dt} \rfloor)

\text{Ans} = \sum_{d=1} d \sum_{t=1} \mu(t) \cdot t^2 \cdot g(\lfloor \frac{n}{dt} \rfloor) \cdot g(\lfloor \frac{m}{dt} \rfloor)

如果此时我们把这个当做最后的式子去进行计算的话，可以发现程序逻辑是：

枚举 $d$ 从 $1$ 到 $n$ 。
在每个 $d$ 内部，枚举 $t$ 从 $1$ 到 $n/d$ 。
计算 $\mu(t) t^2 \dots$

这个调和级数求和的复杂度是

O(N \log N)

。在

N=10^7

的情况下，

10^7 \times \log(10^7) \approx 1.6 \times 10^8

，勉强能跑，但是考虑到我们有大量的取模操作，加上一些神秘的常数，可能导致这道题过不去，最后还需要各种卡常操作，所以我们可以再考虑一下怎么优化最后这个式子。

观察我们上面的式子：

\text{Ans} = \sum_{d=1} d \sum_{t=1} \mu(t) \cdot t^2 \cdot g(\lfloor \frac{n}{dt} \rfloor) \cdot g(\lfloor \frac{m}{dt} \rfloor)

不难发现，只要

d \times t

的积相同，后面的

g

部分就是一样的。那我们为什么不把后面重复算的一堆 提取公因式 合在一起算呢？

我的语言描述能力不太好，可能你听的有点懵，具体得，我们令

T = d \cdot t

。现在我们不再先枚举

d

再枚举

t

，而是直接枚举这个乘积

T

。

现在，我们改写

Ans

的求和顺序：

\text{Ans} = \sum_{T=1}^{\min(n, m)} \left( \sum_{d \cdot t = T} d \cdot \mu(t) \cdot t^2 \right) \cdot g(\lfloor \frac{n}{T} \rfloor) \cdot g(\lfloor \frac{m}{T} \rfloor)

\text{括号里面的} = \sum_{d|T} d \cdot \mu(\frac{T}{d}) \cdot \frac{T^2}{d^2}

这里我们可以消掉一个

d

：

\text{括号里面的} = \sum_{d|T} \frac{T^2}{d} \mu(\frac{T}{d})

我们提取出来一个常数

T

：

\text{括号里面的} = T \sum_{d|T} \frac{T}{d} \mu(\frac{T}{d})

令

k = \frac{T}{d}

，当

d

遍历

T

的所有因子时，

k

也会遍历

T

的所有因子。所以：

\text{括号里面的} = T \sum_{k|T} k \cdot \mu(k)

我们把括号里面的这一堆记进一个函数

F(T)

，原式就变成了:

\text{Ans} = \sum_{T=1}^{\min(n, m)} F(T) \cdot g(\lfloor \frac{n}{T} \rfloor) \cdot g(\lfloor \frac{m}{T} \rfloor)

这样，由于

\lfloor n/T \rfloor

和

\lfloor m/T \rfloor

在一段连续的

T

范围内取值是不变的（整除分块性质），我们可以一次性处理一个区间的

T

，即运用整除分块进行快速计算。

并且，如果我们能预处理出

F(T)

的前缀和

Sum_F(T) = \sum_{i=1}^T F(i)

那么对于一段区间

[l, r]

，它们的贡献就是：

(Sum_F(r) - Sum_F(l-1)) \cdot g(\lfloor n/l \rfloor) \cdot g(\lfloor m/l \rfloor)

现在，我们只剩一个问题：如何快速计算

F(T)

不难发现，

F(T) = T \cdot (Id * \mu)(T)

，其中

*

是狄利克雷卷积，

Id

是单位函数

Id(n)=n

。因为

Id

和

\mu

都是积性函数，它们的卷积也是积性的。

感性理解一下，意思就是：

h(T) = \sum_{k|T} k \mu(k)

是积性的。

所以：

F(T) = T \cdot h(T)

也是积性的。

对于积性函数，我们只需要关心它在质数

p

以及质数的幂

p^k

处的取值，就可以用线性筛

O(n)

扫出来。

这是为什么呢？

根据积性函数的定义：

如果一个函数

f(n)

是积性函数，那么对于任何互质的

a

和

b

（即

\gcd(a, b) = 1

），都有：

f(a \cdot b) = f(a) \cdot f(b)

根据算术基本定理，任何一个大于 1 的正整数

n

都可以唯一分解为：

n = p_1^{a_1} p_2^{a_2} \dots p_k^{a_k}

由于不同的质数幂

p_i^{a_i}

之间一定是互质的，所以：

f(n) = f(p_1^{a_1}) \cdot f(p_2^{a_2}) \cdot \dots \cdot f(p_k^{a_k})

所以只要知道了

f

在所有“质数幂”处的值，就能通过乘法组合出

f

在任何整数处的值。

那么现在考虑

F(T)

在质数幂处的取值：

情况 1： $T$ 是质数 $p$

h(p) = \sum_{k|p} k \mu(k) = 1 \cdot \mu(1) + p \cdot \mu(p) = 1 - p

。

所以：

F(p) = p \cdot (1 - p)

。

情况 2： $T$ 是质数的高次幂 $p^a$ ( $a \ge 2$ ) $h(p^a) = 1 \mu(1) + p \mu(p) + p^2 \mu(p^2) + \dots + p^a \mu(p^a)$ 。

因为当

k \ge 2

时

\mu(p^k) = 0

，所以这个长式子瞬间缩短为：

h(p^a) = 1 - p

。

那么：

F(p^a) = p^a \cdot (1 - p)

。

在代码的线性筛中：

当我们用

i

和质数

p

去筛掉

i \cdot p

时：

情况 A： $i \pmod p \nmid 0$ （即 $i$ 与 $p$ 互质）

根据积性函数的定义，直接相乘：

F(i \cdot p) = F(i) \cdot F(p)

情况 B： $i \pmod p = 0$ （即 $i$ 中已经包含了质因子 $p$ ）

设

i

中包含

p

的最高次数为

p^a

，即

i = M \cdot p^a

（

M

与

p

互质）。

那么

i \cdot p = M \cdot p^{a+1}

。

根据积性：

F(i) = F(M) \cdot F(p^a) = F(M) \cdot p^a (1-p)

F(i \cdot p) = F(M) \cdot F(p^{a+1}) = F(M) \cdot p^{a+1} (1-p)

可以发现，

F(i \cdot p)

刚好就是

F(i)

的

p

倍！

所以：

F(i \cdot p) = F(i) \cdot p

。

我们现在来分析一下优化后的复杂度。

复杂度分析

时间复杂度

线性筛预处理部分，复杂度为 $O(n)$
在线性筛之后，有一个单层循环计算 $F(T)$ 的前缀和,复杂度为 $O(n)$
整除分块计算部分复杂度为 $O(\sqrt{n} + \sqrt{m})$

总体的复杂度还是为

O(n)

\approx 2 \times 10^7

空间复杂度

加上一些零碎的全局变量，总内存占用大约在

210 MB \sim 215 MB

左右。

不过值得注意的是：

根据质数的分布范围分析，实际的 $prime$ 数组不用开这么大（见下面代码）。
如果遇到题目卡常，我们可以运用 $bitset$ $isprime$ ，内存占用会从 $10 MB$ 降到约 $1.25 MB$ 。
$F$ 数组在计算完前缀和之后就不再使用了。所以我们可以直接在 $F$ 数组的基础上原地累加变成 $sum$ 数组，从而省去一个 $long long$ 数组的空间（省下 $80 MB$ ）

但是为了代码的可读性更强，我没有进行以上优化空间复杂度的操作。

以上是这道题的所有推导过程。

参考代码

CPP

#include <bits/stdc++.h>

#define ll long long

using namespace std;

const int mod = 20101009;
const int MAX = 1e7 + 10;

int n, m;

bool isprime[MAX];
int prime[MAX], sub;
ll F[MAX], sum[MAX];

inline void init(int limit)
{
    fill(isprime + 1, isprime + 1 + limit, true);

    F[1] = 1;

    for (int i = 2; i <= limit; ++i)
    {
        if (isprime[i]) { prime[++sub] = i; F[i] = (1ll * i * (1 - i) % mod + mod) % mod; }

        for (int j = 1; j <= sub && i * prime[j] <= limit; ++j)
        {
            isprime[i * prime[j]] = false;

            if (i % prime[j] == 0) { F[i * prime[j]] = F[i] * prime[j] % mod; break; }
            else F[i * prime[j]] = F[i] * F[prime[j]] % mod;
        }
    }

    for (int i = 1; i <= limit; ++i)
    {
        sum[i] = (sum[i - 1] + F[i]) % mod;
    }
}

inline ll g(ll x)
{
    return x * (x + 1) / 2 % mod;
}

int main()
{
    scanf("%d %d", &n, &m);

    init(max(n, m));

    ll ans = 0;

    int limit = min(n, m);

    for (int l = 1, r; l <= limit; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));

        ll res = (sum[r] - sum[l - 1] + mod) % mod;
        ll val = g(n / l) * g(m / l) % mod;

        ans = (ans + res * val % mod) % mod;
    }

    printf("%lld\n", ans);

    return 0;
}

$\textit{\textbf {Eg4} }$

题目描述： P3911 最小公倍数之和

题目简述

求

\sum_{i=1}^N \sum_{j=1}^N \mathrm{lcm}(A_i, A_j)

推导过程

因为

N

很大 (

50000

)，直接

O(N^2)

枚举肯定超时。

但是我们发现数字的大小

A_i

也很小（最大

50000

）。我们可以开一个桶 cnt[x]，表示数字

x

在数组中出现了几次。于是问题转化为枚举数值

i

和

j

（设

V

为最大数值）：

\sum_{i=1}^V \sum_{j=1}^V cnt[i] \cdot cnt[j] \cdot \mathrm{lcm}(i, j)

把

\mathrm{lcm}(i, j)

拆成

\frac{i \cdot j}{\gcd(i, j)}

：

\sum_{i=1}^V \sum_{j=1}^V cnt[i] \cdot cnt[j] \cdot \frac{i \cdot j}{\gcd(i, j)}

像上面那一道题一样，我们先枚举

g = \gcd(i, j)

设

i = g \cdot a, j = g \cdot b

，此时

\gcd(a, b) = 1

。

\sum_{g=1}^V \sum_{a=1}^{V/g} \sum_{b=1}^{V/g} [\gcd(a,b)=1] \cdot cnt[ga] \cdot cnt[gb] \cdot \frac{(ga) \cdot (gb)}{g}

\sum_{g=1}^V g \sum_{a=1}^{V/g} \sum_{b=1}^{V/g} [\gcd(a,b)=1] \cdot (cnt[ga] \cdot a) \cdot (cnt[gb] \cdot b)

这里我们发现

\sum_{a=1}^{V/g} \sum_{b=1}^{V/g} [\gcd(a,b)=1]

很明显可以用莫比乌斯反演进行计算，由于前面这样的计算我们已经进行了很多了，这里我们直接转换

[\gcd(a, b) = 1] = \sum_{d|\gcd(a, b)} \mu(d)

d | \gcd(a, b)

等价于

d|a

且

d|b

。

\sum_{g=1}^V g \sum_{a=1}^{V/g} \sum_{b=1}^{V/g} \left( \sum_{d|a, d|b} \mu(d) \right) \cdot (cnt[ga] \cdot a) \cdot (cnt[gb] \cdot b)

下面又来到了关键的一步，我们同样调换枚举顺序。这里我们考虑先枚举

d

。设

a = d \cdot x

、

b = d \cdot x

。

式子变为

\sum_{g=1}^V g \sum_{d=1}^{V/g} \mu(d) \sum_{x=1}^{\lfloor V/gd \rfloor} \sum_{y=1}^{\lfloor V/gd \rfloor} (cnt[gdx] \cdot dx) \cdot (cnt[gdy] \cdot dy)

这时候我们观察这个式子，可以发现

x

、

y

的枚举是独立的，也就是说，每次循环

x

、

y

的取值是相等的，所以最终的式子可以转换为

Ans = \sum_{g=1}^V g \sum_{d=1}^{V/g} \mu(d) \cdot \left( \sum_{k=1}^{V/(gd)} cnt[k \cdot g \cdot d] \cdot k \cdot d \right)^2

复杂度分析

时间复杂度

运用调和级数进行计算，复杂度大概是

O(V \ln V)

或者

O(V \ln^2 V)

级别

参考代码

CPP

#include <bits/stdc++.h>

#define ll long long

using namespace std;

const int MAX = 50010;

int n, mxaval;
int cnt[MAX];
bool isprime[MAX];
int prime[MAX], sub;
int miu[MAX];

inline void init(int limit)
{
    fill(isprime + 1, isprime + 1 + limit, true);

    miu[1] = 1;

    for (int i = 2; i <= limit; ++i)
    {
        if (isprime[i]) { prime[++sub] = i; miu[i] = -1; }

        for (int j = 1; j <= sub && i * prime[j] <= limit; ++j)
        {
            isprime[i * prime[j]] = false;

            if (i % prime[j] == 0)
            {
                miu[i * prime[j]] = 0;
                break;
            }
            else miu[i * prime[j]] = -miu[i];
        }
    }
}

int main()
{
    scanf("%d", &n);

    mxaval = 0;

    for (int i = 1; i <= n; ++i)
    {
        int x; scanf("%d", &x);
        cnt[x]++;
        mxaval = max(mxaval, x);
    }

    init(mxaval);

    ll ans = 0;

    for (int g = 1; g <= mxaval; ++g)
    {
        ll sumd = 0;

        for (int d = 1; d * g <= mxaval; ++d)
        {
            if (miu[d] == 0) continue;

            ll sumk = 0;
            int T = d * g;

            for (int k = 1; k * T <= mxaval; ++k)
            {
                if (cnt[k * T])
                {
                    sumk += 1ll * cnt[k * T] * k * d;
                }
            }

            if (sumk > 0) sumd += miu[d] * sumk * sumk;
        }

        ans += g * sumd;
    }

    printf("%lld\n", ans);

    return 0;
}

$\textit{\textbf {Eg5} }$

题目描述：P2231 [HNOI2002] 跳蚤

题目简述

题目要求跳蚤通过移动

\sum x_i A_i

（其中

A_i

是卡片上的数字，

x_i

是选了几次这个数字）加上

M

的若干倍（记为

yM

），最终停在

-1

的位置。即方程：

x_1 A_1 + x_2 A_2 + \dots + x_N A_N + y M = -1

要有整数解。

推导过程

根据裴蜀定理（亦或称贝祖定理），上述线性丢番图方程有解的充要条件是：

\gcd(A_1, A_2, \dots, A_N, M) \mid -1

因为 GCD 总是正数，所以条件等价于：

\gcd(A_1, A_2, \dots, A_N, M) = 1

也就是我们要求有多少个序列

(A_1, A_2, \dots, A_N)

，满足

1 \le A_i \le M

，且

\gcd(A_1, \dots, A_N, M) = 1

。

这里我们可以发现这个答案的约束条件很强，直接求解比较难想，又看到了

\gcd

，这时候我们可以稍微往莫比乌斯反演上面想想。

我们设目标函数

f(d)

为满足

\gcd(A_1, \dots, A_N, M) = d

的卡片方案数。我们要求的答案就是

f(1)

。

定义一个辅助函数

F(d)

为满足

d \mid \gcd(A_1, \dots, A_N, M)

的卡片方案数。

由于

1 \le A_i \le M

，这样的数有

M/d

个（即

d, 2d, \dots, \frac{M}{d}d

），而我们可供选择的位置有

N

个，所以：

F(d) = \left( \frac{M}{d} \right)^N

根据莫比乌斯反演公式：

f(n) = \sum_{n|d} \mu(\frac{d}{n}) F(d)

f(1) = \sum_{d|M} \mu(d) F(d) = \sum_{d|M} \mu(d) \left( \frac{M}{d} \right)^N

如果这里我们直接枚举

M

的因子，复杂度会达到

10^8

，并且还要计算

\mu(d)

，很可能会超时。也许你会想到预处理

\mu

？但很明显这是不可取的，会直接

MLE

，所以我们可以考虑怎么稍微优化一下。

观察式子，我们可以发现

\mu(d)

只有在 $d$ 没有平方质因子的时候会产生贡献（不为0），

所以我们可以只关注

M

的质因子。

具体的，我们找出

M

的所有质因子，遍历这些质因子的所有组合，构造因子

d

。然后只需要看

d

里面包含了多少个质因子，算出此时

\mu(d)

的正负，计算相应的贡献即可。

说明

还有一点我想说的，要计算

1 \sim X

之间的质数个数，有一个粗略的计算公式

\pi(X) \approx \frac{X}{\ln X}

在计算时，我们可以粗略的认为

\ln X \approx 2.3 \log_{10} X

。

如果想要更深层的了解

\pi(X)

,可以参考

1.《初等数论》（潘承洞、潘承彪著，北京大学出版社，第四版）、

2.《初等数论》（陈景润著，哈尔滨工业大学出版社）

3.《数论：概念和问题》（蒂图·安德雷斯库著，哈尔滨工业大学出版社）等著作，里面都有关于此不错的讲解。

此外，我们也可以通过计算质数的乘积来粗略计算。

例如对于这道题，

2 \times 3 \times 5 \times 7 \times 11 \times 13 \times 17 \times 19 \times 23 = 223,092,870

，只需要9个从小到大的质数，就已经超过了

10^8

，所以我们只需要开大小为

10

的数组就好了。

参考代码

CPP

#include <bits/stdc++.h>

#define ll long long

using namespace std;

const int MAX = 50;

int n, m;

int prime[MAX], sub;

ll ans = 0;

inline ll qmi(ll a, int b)
{
    ll res = 1;

    while (b)
    {
        if (b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }

    return res;
}

void dfs(int idx, int d, int cnt)
{
    if (idx > sub)
    {
        ll term = qmi((ll)m / d, n);
        ans += (cnt % 2 == 0 ? term : -term);
        return;
    }

    dfs(idx + 1, d, cnt); dfs(idx + 1, d * prime[idx], cnt + 1);
}

int main()
{
    scanf("%d %d", &n, &m);

    int temp = m;

    for (int i = 2; i * i <= temp; ++i)
    {
        if (temp % i == 0)
        {
            prime[++sub] = i;
            while (temp % i == 0) temp /= i;
        }
    }

    if (temp > 1) prime[++sub] = temp;

    dfs(1, 1, 0);

    printf("%lld\n", ans);

    return 0;
}

$\textit{\textbf {Eg6} }$

题目描述 P3312 [SDOI2014] 数表

题目简述

求

Ans

对

2^{13}

取模的结果。其中，

\text{Ans} = \sum_{i=1}^n \sum_{j=1}^m [f(\gcd(i, j)) \le a] \cdot f(\gcd(i, j))

推导过程

忽略

f(\gcd(i, j)) \le a

的限制，先推导通项公式。

令

d = \gcd(i, j)

，枚举

d

：

\sum_{d=1}^{\min(n, m)} f(d) \sum_{i=1}^n \sum_{j=1}^m [\gcd(i, j) = d]

\sum_{d=1}^{\min(n, m)} f(d) \sum_{i=1}^{\lfloor n/d \rfloor} \sum_{j=1}^{\lfloor m/d \rfloor} [\gcd(i, j) = 1]

利用莫比乌斯反演，

\sum_{d=1}^{\min(n, m)} f(d) \sum_{k=1}^{\min(\lfloor n/d \rfloor, \lfloor m/d \rfloor)} \mu(k) \lfloor \frac{n}{dk} \rfloor \lfloor \frac{m}{dk} \rfloor

令

T = d \cdot k

，则

d

是

T

的约数。我们要枚举

T

：

\sum_{T=1}^{\min(n, m)} \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \sum_{d|T} f(d) \cdot \mu(\frac{T}{d})

现在考虑限制条件

f(d) \le a

。

公式变为：

\text{Ans} = \sum_{T=1}^{\min(n, m)} \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \left( \sum_{d|T, f(d) \le a} f(d) \cdot \mu(\frac{T}{d}) \right)

内层的求和

\sum_{d|T, f(d) \le a} f(d) \cdot \mu(\frac{T}{d})

依赖于

a

。由于有多组查询，且

a

各不相同，我们不能针对每个

a

重新计算。

我们可以考虑用离线 + 树状数组解决

排序：将所有查询按 $a$ 从小到大排序。同时，将所有数字 $1 \dots 10^5$ 按照其约数和 $f(d)$ 从小到大排序。
动态更新：遍历排序后的查询。对于当前的查询限制 $a$ ，将所有满足 $f(d) \le a$ 的数字 $d$ 加入到数据结构中。

当我们加入一个数字 $d$ 时，它会对所有 $T$ (其中 $T$ 是 $d$ 的倍数) 产生贡献。
贡献值为 $f(d) \cdot \mu(T/d)$ 。
我们需要一个支持单点修改、区间求和的数据结构来维护 $g(T) = \sum f(d) \mu(T/d)$ 的前缀和。毫无疑问，树状数组是最佳选择。

整除分块：对于每个查询，利用维护好的树状数组，结合 $\lfloor n/T \rfloor$ 和 $\lfloor m/T \rfloor$ 的整除分块（数论分块）在 $O(\sqrt{n} \log n)$ 时间内算出答案。

复杂度

时间复杂度

线性筛预处理 $\mu$ 和 $f$ (即 $\sigma_1$ )： $O(N)$ 。
排序： $O(N \log N + Q \log Q)$ 。
离线处理 + 树状数组更新：每个数 $d$ 会更新其所有倍数，总复杂度 $\sum \frac{N}{d} = O(N \log N)$ 。
查询：每个查询 $O(\sqrt{N} \cdot \log N)$ 。
总复杂度： $O(N \log N + Q\sqrt{N} \log N)$ ，通过这道题绰绰有余

参考代码

CPP

#include <bits/stdc++.h>

#define ll long long

using namespace std;

const ll MOD = 1 << 31;
const int MAX = 1e5 + 10;

int Q;

int prime[MAX], sub;
bool isprime[MAX];

int miu[MAX];
int sigma[MAX]; // 因子和
int low[MAX]; // i 的最小质因子p的最大幂次的等比数列和
int smallp[MAX]; // 记录最小质因子
int bit[MAX]; // 树状数组

struct node
{
    int id, val;
}f[MAX];

struct ASK
{
    int n, m, a, id;
}q[MAX];

int ans[MAX];

inline bool compare1(node x, node y)
{
    return x.val < y.val;
}

inline bool compare2(ASK x, ASK y)
{
    return x.a < y.a;
}

inline void init(int limit)
{
    fill(isprime + 1, isprime + 1 + limit, true);

    miu[1] = 1; sigma[1] = 1;

    for (int i = 2; i <= limit; ++i)
    {
        if (isprime[i])
        {
            prime[++sub] = i;
            miu[i] = -1;
            sigma[i] = i + 1; // 质数的因数只有1和他自己
            low[i] = i + 1;
            smallp[i] = i;
        }

        for (int j = 1; j <= sub && i * prime[j] <= limit; ++j)
        {
            isprime[i * prime[j]] = false;

            if (i % prime[j] == 0)
            {
                miu[i * prime[j]] = 0;
                smallp[i * prime[j]] = smallp[i] * prime[j];
                low[i * prime[j]] = low[i] * prime[j] + 1;
                sigma[i * prime[j]] = sigma[i / smallp[i]] * low[i * prime[j]];
                break;
            }
            else
            {
                miu[i * prime[j]] = -miu[i];
                smallp[i * prime[j]] = prime[j];
                low[i * prime[j]] = prime[j] + 1;
                sigma[i * prime[j]] = sigma[i] * sigma[prime[j]];
            }
        }
    }

    for (int i = 1; i <= limit; ++i)  f[i].id = i, f[i].val = sigma[i];
}

inline int lowbit(int x)
{
    return x & -x;
}

inline void add(int x, int val)
{
    for (; x < MAX; x += lowbit(x)) bit[x] += val;
}

inline int query(int x)
{
    int res = 0;

    for (; x > 0; x -= lowbit(x)) res += bit[x];

    return res;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0); cout.tie(0);

    cin >> Q;

    init(100000);

    for (int i = 1; i <= Q; ++i)
    {
        cin >> q[i].n >> q[i].m >> q[i].a;
        q[i].id = i;
    }

    sort(f + 1, f + 100001, compare1);
    sort(q + 1, q + 1 + Q, compare2);

    int pos = 1;

    for (int i = 1; i <= Q; ++i)
    {
        while (pos <= 100000 && f[pos].val <= q[i].a)
        {
            int d = f[pos].id;
            for (int T = d; T <= 100000; T += d)
            {
                add(T, f[pos].val * miu[T / d]);
            }

            pos++;
        }

        int n = q[i].n, m = q[i].m;

        if (n > m) swap(n, m);

        ll res = 0;

        for (int l = 1, r; l <= n; l = r + 1)
        {
            r = min(n / (n / l), m / (m / l));
            res += (n / l) * (m / l) % MOD * (query(r) - query(l - 1)) % MOD;
        }

        ans[q[i].id] = res % MOD;
    }

    for (int i = 1; i <= Q; ++i) cout << ans[i] << endl;

    return 0;
}

说明

代码中的

low[MAX]

存的是数字

i

的最小质因子

p

的等比数列之和。

简单来说，如果一个数

i

的质因数分解为：

i = p_1^{c_1} \cdot p_2^{c_2} \cdots p_k^{c_k}

假设

p_1

是最小的质因子，那么 low[i] 存的就是：

\text{low}[i] = 1 + p_1 + p_1^2 + \dots + p_1^{c_1}

可以利用这个算数和函数 $\sigma(i)$ 。

结语

以上就是 莫比乌斯反演 的全部讲解啦，到这里你应该对莫比乌斯反演这个神奇的工具有了一定的了解了。最后给大家推荐一个自己建的题单：莫比乌斯反演的运用。里面有一些关于莫比乌斯反演的好题，还是值得一做的（保证可以用莫比乌斯反演完成）。

感谢大家耐心看完，希望这篇文章对你有所帮助。

完结撒花 QWQ

文章操作

前置知识

莫比乌斯函数

定义

性质

线性筛代码