「算法笔记」莫比乌斯反演

莫比乌斯反演是数论中的一种小技巧，常用于解决一些数论函数的求和问题。通常做法是使用 $\mu$ 函数来将 $[n=1]$ 转换成 $\sum_{d|n}\mu(i)$ 的形式。

推荐博客：algocode 算法博客

$\mu$ 的定义：

• 如果 $n$ 有至少一个平方因子，那么 $\mu(n)=0$
• 如果 $n$ 没有一个平方因子，且有 $2k+1$ 个素因子，那么 $\mu(n)=-1$
• 如果 $n$ 没有一个平方因子，且有 $2k$ 个素因子，那么 $\mu(n)=1$

$\mu$ 的性质：

• $\mu$ 是积性函数。
• 根据这一性质，我们可以使用线性筛在 $\Theta(n)$ 的时间中筛出 $i\in[1,n]$ 中的所有 $\mu(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];
    }
}

• $\sum_{d|n}\mu(i)=1$
• 证明见 algocode 算法博客。
• 利用这一性质，我们就可以解决一些简单的问题了。

例题一：BZOJ 2301 Problem B

给定 $a,b,c,d,e$，求 $\sum_{i=a}^{b}\sum_{j=c}^{d}[gcd(i,j)=k]$ 的值。
$1\le a\le b,c\le d,k\le 50000$，$50000$ 组询问。

可以将询问拆成四个前缀询问。
问题变成了：求 $\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)=k]$ 的值。

将 $i,j$ 都除以 $k$，上式等于 $\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{k}\rfloor}[gcd(i,j)==1]$

利用 $\mu$ 函数的第二个性质，上式等于 $\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{d|i,d|j}\mu(d)$

将枚举 $d$ 移到前面，并令 $i=i'd,j=j'd$，
则上式等于 $\sum_{d=1}^{n}\mu(d)\sum_{i'=1}^{\lfloor\frac{n}{kd}\rfloor}\sum_{j'=1}^{\lfloor\frac{n}{kd}\rfloor}1=\sum_{d=1}^{n}\mu(d)\lfloor\frac{n}{kd}\rfloor^2$

而 $\lfloor\dfrac{n}{i}\rfloor$ 只有 $\Theta(\sqrt n)$ 种取值。直接预处理 $\mu$ 函数的前缀和，然后数论分块即可。

#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;
}

例题二：BZOJ 2154 Crash 的数字表格

求：$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)$。$n,m\le 10^7$

发现 $lcm(i,j)=\dfrac{ij}{gcd(i,j)}$。

考虑枚举 $gcd(i,j)=d$。问题就变成了：

$$\sum_{d=1}^{min(n,m)}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=1]\cdot ij$$

令 $f(n,m)$ 表示 $\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=1]\cdot ij$，我们发现可以使用莫比乌斯反演得出原式等于：

$$\sum_{d=1}^{min(n,m)}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij$$

$$=\sum_{d=1}^{min(n,m)}\mu(d)\dfrac{\lfloor\frac{n}{d}\rfloor(\lfloor\frac{n}{d}\rfloor+1)}{2}\dfrac{\lfloor\frac{m}{d}\rfloor(\lfloor\frac{m}{d}\rfloor+1)}{2}$$
这个式子可以使用数论分块来做。

然后，回到原问题。我们再使用一次数论分块，这样就可以在 $\Theta((\sqrt n)^2)=\Theta(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;
}