【Coel.学习笔记】莫比乌斯反演

闲话

记得在差不多一年前写扩展欧拉定理的时候我提了一句

这周终于把古代猪文搞定了，数论这块的内容就只剩个博弈论了
别提莫比乌斯反演之类的东西，我不想搞

甚至刚开始写的时候还笔误把莫反写成了莫队……

转眼一年过去了，来填上这个大坑吧！

前言

莫比乌斯反演是一个基于莫比乌斯函数 \(\mu(n)\) 的反演算法。对于某些求和函数，如果直接求解的时间复杂度很高，而倍数或因数之和求解较为容易，则可以通过莫比乌斯反演简化运算，提高效率。

前置知识

这些前置知识可以算到莫比乌斯反演的一部分，不过也可以作为一个独立的内容（在进阶指南里其实就是这样划分的）。

数论分块

差不多相当于一个定理，这里就不多提了。建议自行查资料然后做做题，熟练后再看莫反。

莫比乌斯函数的定义

记正整数 \(x\) 的质因数分解结果为 \(x=p_1^{\alpha_1}p_2^{\alpha_2}...p_k^{\alpha_n}\)（\(p_i\) 为质数且 \(\alpha_i\geq 1\)），则

\[\mu(x)=\begin{cases} 0 & \exists i,\alpha_i>1\text{（即含有平方因子）} \\ (-1)^k & \forall i,\alpha_i=1\text{ （即不含有平方因子） } \end{cases} \]

特殊地，\(\mu(1)=1\)。

莫比乌斯函数是一个积性函数，这意味着我们可以用埃氏筛或者线性筛预处理。

莫比乌斯函数的性质

记 \(n\) 的约数的莫比乌斯函数之和 \(S_n = \sum_{d | n}\mu(d)\)，则

\[S_n=\begin{cases}1 & n=1 \\ 0 & \text{otherwise.}\end{cases} \]

\(n=1\) 显然正确，其他情况证明利用二项式定理容易得到，这里略去。

顺带一提，这个性质可以说明 \(\mu\) 在狄利克雷生成函数中为黎曼函数 \(\zeta\) 的逆元，不过太高深没啥用（

反演定理

若对于函数 \(F(n),f(n)\)，满足 \(F(n)=\sum_{d|n}f(d)\)，则 $$f(n)=\sum_{d|n}\mu(d)F(\dfrac{n}{d})$$

下面是证明。我们证明得到的式子左边等于右边即可：

\[\begin{aligned}\sum_{d|n}\mu(d)F(\dfrac{n}{d}) &= \sum_{d|n}\mu(d)\sum_{i|\frac{n}{d}}f(i) & \text{【代入函数 F】}\\ &=\sum_{i|n}f(i)\sum_{d|\frac{n}{i}}\mu(d)& \text{【进行数论变换】}\\ &=f(n)& \text{【利用前面提到的性质】} \end{aligned} \]

\(\text{Q.E.D.}\)

---

另一种反演定理的形式为：若 \(F(n)=\sum_{n|d}f(d)\)，则 \(f(n)=\sum_{n|d}\mu(\dfrac{d}{n})F(d)\)。

类比上一个的证明过程也不难得证：

\[\begin{aligned}\sum_{n|d}\mu(\dfrac{d}{n})F(d) &=\sum_{n|d}\mu(\dfrac{d}{n})\sum_{d|i}f(i) \\ &=\sum_{n|i}f(i)\sum_{d^\prime|\frac{i}{n}}\mu(d^\prime)\\ &= \sum_{n|i}f(i)S(\dfrac{i}{n})\\ &=f(n)\end{aligned} \]

\(\text{Q.E.D.}\)

利用这两个定理，我们可以把从 \(f(n)\) 到 \(F(n)\) 翻转过来求解，这个过程就叫做莫比乌斯反演。

例题

[POI2007]ZAP-Queries

洛谷传送门
求解

\[\sum^a_{i=1}\sum^b_{j=1}[gcd(i,j)=d] \]

多组数据。

解析：这道题的推导过程比较简单，所以放在最前面。

按照莫比乌斯反演的定理，我们需要定义两个函数 \(f(k),F(k)\)。对于这道题来说：

\[f(k)=\sum^a_{i=1}\sum^b_{j=1}[gcd(i,j)=k] \]

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

则根据反演定理可以直接得到

\[f(k)=\sum_{k|d}\mu(\dfrac{d}{k})F(k) \]

改变枚举顺序

\[f(k)=\sum_{i=1}^{\min(a,b)}\mu(i)\lfloor\dfrac{a}{ik}\rfloor\lfloor\dfrac{b}{ik}\rfloor \]

数论分块处理即可。

#include <iostream>
#include <cstring>
#include <bitset>

using namespace std;

const int maxn = 5e4 + 10;

namespace __Coel_FastIO {
    const int SIZE = 1 << 20;
    static char buf_r[SIZE], *p1(buf_r), *p2(buf_r);
    static char buf_w[SIZE], *p3(buf_w), *p4 = (buf_w + SIZE);
    static std::streambuf *inbuf = std::cin.rdbuf(), *outbuf = std::cout.rdbuf();
    inline static void gc(char &c) {
        c = p1 == p2 && (p2 = (p1 = buf_r) + inbuf->sgetn(buf_r, SIZE), p1 == p2) ? EOF : *p1++;
    }
    inline static void pc(char c) {
        p3 == p4 ? outbuf->sputn(p3 = buf_w, SIZE), *p3++ = c : *p3++ = c;
    }

    inline static int Flush() {
        return outbuf->sputn(buf_w, p3 - buf_w), 0;
    }

    struct IOdesu {
        template<typename T> IOdesu& operator >>(T &x) {
            x = 0;
            int f = 1;
            char ch;
            gc(ch);
            while (ch < '0' || ch > '9') {
                if (ch == '-') f = -1;
                gc(ch);
            }
            while (ch >= '0' && ch <= '9') x = x * 10 + ch - '0', gc(ch);
            return x *= f, *this;
        }

        IOdesu& operator <<(const char ch) {
            return pc(ch), *this;
        }

        template<typename T> IOdesu& operator <<(T x)  {
            if (x < 0) x = -x, pc('-');
            static int buf[35];
            int top = 0;
            do {
                buf[top++] = x % 10;
                x /= 10;
            } while (x);
            while (top) pc(buf[--top] + '0');
            return *this;
        }

    } qwq;
}

using namespace __Coel_FastIO;
using namespace std;

int T, cnt;
int mu[maxn], s[maxn], prime[maxn];
bitset<maxn> vis;

inline int get(int k, int x) {
    return k / (k / x);
}

void initMu() {
    mu[1] = 1;
    for (int i = 2; i < maxn; i++) {
        if (!vis[i]) prime[cnt++] = i, mu[i] = -1;
        for (int j = 0; i * prime[j] < maxn; j++) {
            vis[prime[j] * i] = true;
            if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (int i = 1; i < maxn; i++) s[i] = s[i - 1] + mu[i];
}

int main(void) {
    initMu();
    qwq >> T;
    while (T--) {
        int a, b, d, k, res = 0;
        qwq >> a >> b >> d;
        k = min(a, b);
        for (int l = 1, r; l <= k; l = r + 1) {
            r = min(get(a, l), get(b, l));
            res += (a / (l * d)) * (b / (l * d)) * (s[r] - s[l - 1]);
        }
        qwq << res << '\n';
    }
    return Flush();
}

[HAOI2011] problem b

洛谷传送门
求解

\[\sum^b_{i=a}\sum^d_{j=c}\gcd (i,j) \]

多组数据。
解析：这道题是我之前出过的一道题的简易版（可恶，居然是重题），推导过程参考那篇题解就行了，这里不再赘述。

[SDOI2015] 约数个数和

洛谷传送门
求解

\[\sum^n_{i=1}\sum^m_{j=1}\text{d} (i\times j) \]

同样多组数据，其中 \(\text{d}(i)\) 表示 \(i\) 的约数个数。

解析：这道题直接对着 \(\text{d}\) 求很难，考虑利用 \(i\times j\) 的性质。

\[\text{Lemma:} \text{d}(i\times j)=\sum_{x|i}\sum_{y|j} [\gcd (x,y)=1] \]

这个引理的证明可以利用分解质因数得到，这里略过。

根据这个原理，原式转化为

\[\sum^n_{i=1}\sum^m_{j=1}\sum_{x|i}\sum_{y|j} [\gcd (x,y)=1] \]

改变求和顺序有

\[\sum^n_{x=1}\sum^m_{y=1}\lfloor\dfrac{n}{x}\rfloor\lfloor\dfrac{m}{y}\rfloor[\gcd(x,y)=1] \]

接下来进行反演。令 $$F(k)=\sum^n_{i=1} \sum^m_{j=1}\lfloor\dfrac{n}{i}\rfloor\lfloor\dfrac{m}{j}\rfloor[\gcd(x,y)=k]$$

\[f(k)=\sum^n_{i=1} \sum^m_{j=1}\lfloor\dfrac{n}{i}\rfloor\lfloor\dfrac{m}{j}\rfloor[k|\gcd(i,j)] \]

那么有

\[f(k)=\sum_{k|d}f(d) \]

对 \(F(k)\) 提取 \(k\) 可以得到

\[F(k)=\sum^{n/k}_{i=1}\sum^{m/k}_{j=1}\lfloor\dfrac{n}{ik}\rfloor\lfloor\dfrac{m}{jk}\rfloor \]

根据之前定义可知答案为 \(F(1)\)。又由上面得到的反演定理可知

\[F(k)=\sum_{k|d}\mu(\dfrac{d}{k})f(d) \]

把 \(k=1\) 代入可得

\[F(1)=\sum^{\min(n,m)}_{i=1}\mu(i)f(i) \]

为了便于计算，定义函数

\[h(k)=\sum_{i=1}^{k}\lfloor{\dfrac{k}{i}}\rfloor \]

那么把 \(h(k)\) 预处理出来，就可以 \(O(1)\) 求解 \(f(i)\)，再把 \(\mu(i)\) 的前缀和求出来，剩下的交给数论分块即可。时间复杂度为 \(O(n\sqrt n +T\sqrt n)\)，前者为数论分块预处理 \(h(k)\) 的时间复杂度。

（顺带一提，这代码在 loj 上跑得飞快，洛谷上 \(999\text{ms}\)，acwing 上必须开 O3 优化才能过……）

void initMu() {
    mu[1] = 1;
    for (int i = 2; i < maxn; i++) {
        if (!vis[i]) prime[cnt++] = i, mu[i] = -1;
        for (int j = 0; i * prime[j] < maxn; j++) {
            vis[prime[j] * i] = true;
            if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (int i = 1; i < maxn; i++) s[i] = s[i - 1] + mu[i];
    for (int i = 1; i < maxn; i++) {
        for (int l = 1, r; l <= i; l = r + 1) {
            r = min(i, get(i, l));
            h[i] += (r - l + 1) * (i / l);
        }
    }
}

signed main(void) {
    initMu();
    qwq >> T;
    while (T--) {
        int n, m, res = 0, k;
        qwq >> n >> m;
        k = min(n, m);
        for (int l = 1, r; l <= k; l = r + 1) {
            r = min(k, min(get(n, l), get(m, l)));
            res += (s[r] - s[l - 1]) * h[n / l] * h[m / l];
        }
        qwq << res << '\n';
    }
    return Flush();
}