莫比乌斯反演

发布时间 2023-08-27 07:16:38作者: -白简-

前置知识

狄利克雷卷积和积性函数

数论分块学习笔记

和式简单推倒技巧

莫比乌斯函数

定义

\[\mu(n)= \begin{cases} 1&n=1\\ 0&n\text{ 含有平方因子}\\ (-1)^k&k\text{ 为 }n\text{ 的本质不同质因子个数}\\ \end{cases} \]

设 \(n=\prod^{k}_{i=1}p_i^{c_i}\),其中 \(p_i\) 为 \(n\) 的质因子,\(c_i \geq 1\)。

性质

莫比乌斯函数是积性函数,有以下性质

\[\sum_{d \mid n}\mu(d)=\left[n=1 \right] \]

常用于莫比乌斯反演中。

证明

  1. 当 \(n=1\) 时显然成立;
  2. 当 \(n \neq 1\) 时,设

\[n=\prod^{k}_{i=1}p_i^{c_i} \]

那么有

\[\sum_{d \mid n}\mu(d)=\sum^{n}_{i=0}\dbinom{k}{i}\cdot(-1)^i \]

考虑我们选 \(i\) 个本质不同的质因子,有 \(\binom{k}{i}\) 种方法,其中每个数的莫比乌斯函数值为 \((-1)^i\),所以这样计算出原式的贡献。

那么我们需要证明

\[\sum^{n}_{i=0}\dbinom{k}{i}\cdot(-1)^i=0 \]

考虑二项式定理

\[\sum^{n}_{i=0}\dbinom{k}{i}\cdot(-1)^i=\sum^{n}_{i=0}\dbinom{n}{i} \cdot (-1)^i \cdot 1^{n-i}=(-1+1)^n=0 \]

证毕。

莫比乌斯反演

反演公式

设 \(f(n),g(n)\) 是定义在正整数集上的两个函数。

形式 \(1\)

若有

\[f(n)=\sum_{d \mid n}g(d)=\sum _{d \mid n}g(\dfrac{n}{d})\]

则有

\[g(n)=\sum_{d \mid n}f(\dfrac{n}{d})\mu(d)=\sum_{d \mid n}f(d)\mu(\dfrac{n}{d}) \]

证明

考虑对下式进行变换

\[\sum_{d \mid n}\mu(d)f(\dfrac{n}{d}) \]

代入定义,有

\[\sum_{d \mid n}\mu(d)f(\dfrac{n}{d})=\sum_{d\mid n}\mu(d)\sum_{t\mid \frac{n}{d}}g(t) \]

可以交换求和号

\[\sum_{t\mid n}g(t)\sum_{d\mid \frac{n}{t}}\mu(d) \]

详细解释 对于上面这一步变换,有以下考虑。

由于 \(d \mid n\),存在 \(t\) 满足 \(t \mid \!\frac{n}{d}\)。那么对于每一个符合条件的 \(t\) 都有与之对应的 \(d\),那么就可以转化为

\[\sum_{dt \mid n}\mu(d)g(t) \]

考虑原式,相当于先枚举符合条件的 \(d\),再枚举在此时 \(d\) 基础上符合条件的 \(t\),我们发现,换成现在的式子后,对答案的统计仍然不重不漏。

还可以换一种方式考虑,对于原式,我们可以转化为

\[\sum_d \mu(d) \left[ d \mid n\right]\sum_t g(t)\left[dt \mid n\right] \]

显然,对于两个限制条件,我们可以直接消去第一个条件,再将艾弗森括号转换进和式,得到

\[\sum_{dt \mid n}\mu(d)g(t) \]

得到这个式子,我们就可以直接交换 \(\mu(d)\) 和 \(g(t)\) 的位置了,并可以将和式中的一个条件拆开成为两个和式

\[\sum_{t \mid n}g(t)\sum_{d \mid \frac{n}{t}}\mu(d) \]

又由 \(\sum_{d\mid n}\mu(d)=\left[ n=1 \right]\) 可以得到

\[=\sum_{t \mid n}g(t)\left[ \frac{n}{t}=1 \right] \]

把艾弗森括号里的部分转化一下

\[=\sum_{t \mid n}g(t)\left[ n=t \right] \]

只有当 \(n=t\) 时,该式子才有值,即

\[=g(n) \]

证毕。

形式 \(2\)

若有

\[f(n)=\sum_{n \mid d}g(d)=\sum_{n \mid d}g(\frac{n}{d}) \]

则有

\[g(n)=\sum_{n \mid d}\mu(\frac{d}{n})f(d)=\sum_{n \mid d}\mu(d)f(\frac{d}{n}) \]

证明与形式 \(1\) 的证明类似,不再展开说了。

例题

POI2007 ZAP-Queries

题目大意

给定 \(n\) 和 \(m\),求

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

证明

为了方便处理,我们钦定 \(n \leq m\)。

首先考虑 \(d\) 同时是 \(i\) 和 \(j\) 的最大公约数,所以满足 \(d \ | \ i\) 且 \(d \ | \ j\)。

于是我们可以把 \(d\) 提出来。

\[\sum_{i = 1}^{\left\lfloor\\\frac{n}{d}\right\rfloor}\sum_{j = 1}^{\left\lfloor\\\frac{m}{d}\right\rfloor} \left[ \gcd(i,j) = 1\right] \]

然后我们可以利用莫比乌斯函数的结论,变为

\[\sum_{i = 1}^{\left\lfloor\\\frac{n}{d}\right\rfloor}\sum_{j = 1}^{\left\lfloor\\\frac{m}{d}\right\rfloor}\sum_{t |\gcd(i,j)} \mu(t) \]

此时,实际上 \(t\) 是随意枚举的,只要保证 \(t\ | \gcd(i,j)\) 就可以。既然 \(t\) 不受到前面两个和式的限制,就可以向前提,先枚举 \(t\),又考虑到满足 \(t\ | \gcd(i,j)\) 的最大的值也就是 \(\min(\left\lfloor\\\frac{n}{d}\right\rfloor,\left\lfloor\\\frac{m}{d}\right\rfloor)\),所以 \(t\) 就枚举到这就行。

又考虑对 \(\mu(t)\) 产生限制的只有 \(t\),那么就可以把她也向前提前面去,放在 \(t\) 的和式后面。我们也可以顺便把其他和式的限制条件放在她后面。

\[\sum_{t = 1}^{\left\lfloor\\\frac{n}{d}\right\rfloor}\mu(t) \sum_{i = 1}^{\left\lfloor\frac{n}{d} \right\rfloor}\left[t \ | \ i \right]\sum_{j=1}^{\left\lfloor\\\frac{m}{d}\right\rfloor} \left[t \ | \ j \right] \]

然后我们可以发现,观察后面两个和式。例如第二个和式,满足 \(\left[t \ | \ i \right]\) 的 \(i\) 最多有 \(\left\lfloor\\\frac{n}{dt}\right\rfloor\) 个。

那么我们这个式子又可以改变,变成下面这个样子。

\[\sum_{t = 1}^{\left\lfloor\\\frac{n}{d}\right\rfloor}\mu(t) \left\lfloor\\\frac{n}{dt}\right\rfloor \left\lfloor\\\frac{m}{dt}\right\rfloor \]

对于后面两个向下取整的部分,我们可以整除分块求,前面的部分我们可以预处理求出莫比乌斯函数。

Code
#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cmath>

using namespace std;

#define int long long

const int N = 50050;

int T,n,m,d;

int pri[N],p[N],cnt;
int miu[N],mu[N];

void Ready() {
    p[1] = 1;
    mu[1] = 1;

    for(int i = 2;i < N; i++) {
        if(p[i] == 0) {
            pri[++cnt] = i;
            mu[i] = -1;
        }

        for(int j = 1;j <= cnt && i * pri[j] < N; j++) {
            p[i * pri[j]] = 1;

            if(i % pri[j] == 0)
                break;
            
            mu[i * pri[j]] = -mu[i];
        }
    }

    for(int i = 1;i < N; i++) 
        miu[i] = mu[i] + miu[i - 1];
    return ;
}

int ans = 0;

signed main() {
    cin >> T;

    Ready();

    while(T--) {
        cin >> n >> m >> d;
        if(n > m)
            swap(n,m);
        n /= d;
        m /= d;
        int l = 1,r = 0;
        while(l <= n) {
            r = min(n / (n / l),m / (m / l));
            ans += (n / l) * (m / l) * (miu[r] - miu[l - 1]);
            l = r + 1;
        }

        cout << ans << "\n";
        ans = 0;
    }
    return 0;
}

YY的GCD

题目大意

给定 \(n\) 和 \(m\),求

\[\sum_{i=1}^n\sum^{m}_{j=1}\left[\gcd(i,j) \in \mathbb{P} \right] \]

即求最大公约数为质数的数对个数。

思路

为了方便处理,我们还是钦定 \(n \leq m\)。

我们可以枚举一个 \(d\),把原式转化为

\[\sum_{d=1}^{n} \left[ d \in \mathbb{P} \right]\sum_{i=1}^n\sum^{m}_{j=1}\left[\gcd(i,j) = d \right] \]

利用莫比乌斯函数性质,有

\[\sum_{d=1}^{n} \left[ d \in \mathbb{P} \right]\sum_{i=1}^{\left\lfloor \frac{n}{d}\right\rfloor} \sum^{\left\lfloor \frac{m}{d}\right\rfloor}_{j=1}\sum^{}_{t \mid \gcd(i,j)}\mu(t) \]

交换求和号,令 \(i=\frac{i}{t},j=\frac{j}{t}\),有

\[\begin{aligned} &= \sum_{d=1}^{n} \left[ d \in \mathbb{P} \right] \sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{t=1}\mu(t) \sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{i=1}\left[i \mid t\right] \sum^{\left\lfloor \frac{m}{d}\right\rfloor}_{j=1}\left[j \mid t\right] \\ &= \sum_{d=1}^{n} \left[ d \in \mathbb{P} \right] \sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{t=1}\mu(t) \sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{i=1}\left[i \mid t\right] \sum^{\left\lfloor \frac{m}{d}\right\rfloor}_{j=1}\left[j \mid t\right] \\ &= \sum_{d=1}^{n} \left[ d \in \mathbb{P} \right] \sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{t=1} \mu(t) {\left\lfloor \frac{n}{dt}\right\rfloor} {\left\lfloor \frac{m}{dt}\right\rfloor} \\ \end{aligned} \]

式子推到这里,如果我们去枚举 \(d\),肯定会 TLE,我们令 \(T = dt\),有

\[\sum^{n}_{T=1}{\left\lfloor \frac{n}{T}\right\rfloor}{\left\lfloor \frac{m}{T}\right\rfloor}\sum^{}_{d \mid T}\left[ d \in \mathbb{P} \right]\mu(\frac{T}{d}) \]

我们发现,后面那部分可以预处理。

考虑每个质数 \(d\),对于它的倍数 \(T\),其值加上 \(\mu(\frac{T}{d})\)。

Code
#include <bits/stdc++.h>

using namespace std;

const int N = 10050000;

int T,n,m;

bool p[N];

int mu[N],sum[N],f[N];

int prime[N],cnt;

void Seive() {
    mu[1] = 1;

    for(int i = 2;i <= 10000000; i++) {
        if(!p[i]) {
            cnt ++;
            prime[cnt] = i;
            mu[i] = -1;
        }

        for(int j = 1;i * prime[j] <= 10000000 && j <= cnt; j++) {
            p[i * prime[j]] = 1;

            if(i % prime[j] == 0)
                break;
            
            mu[i * prime[j]] = - mu[i];
        }
    }
    
    for(int i = 1;i <= cnt; i++) 
        for(int j = 1;prime[i] * j <= 10000000; j++) 
            f[prime[i] * j] += mu[j];

    for(int i = 1;i <= 10000000; i++)
        sum[i] = sum[i - 1] + f[i];
}

long long H(int a,int b) {
    long long res = 0;

    if(a == 0)
        return 0;

    int l = 1,r = 0;
    while(l <= n) {
        r = min(a / (a / l),b / (b / l));
        res += (sum[r] - sum[l - 1]) * (long long)(a / l) * (long long)(b / l);
        l = r + 1;
    }
    return res;
}

int main() {
    ios::sync_with_stdio(false);

    cin >> T;

    Seive();
    
    while(T --> 0) {
        cin >> n >> m;

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

        cout << H(n,m) << "\n";
    }
    return 0;
}

SDOI2014 数表

题目大意

有一张 \(n\times m\) 的数表,其第 \(i\) 行第 \(j\) 列(\(1\le i\le n\),\(1\le j\le m\))的数值为能同时整除 \(i\) 和 \(j\) 的所有自然数之和。给定 \(a\),计算数表中不大于 \(a\) 的数之和。

(\(T \leq 2 \times 10^4,1 \leq n,m \leq 10^5,a\leq 10^9\))。

思路

形式化题意:

给定 \(n,m,a\),求

\[\sum^{n}_{i=1} \sum^{m}_{j=1} \sigma_1(\gcd\left(i, j\right)) \left[ \sigma_1(\gcd\left(i, j\right)) \le a\right] \]


为了方便处理,我们还是钦定 \(n \leq m\)。

先不考虑 \(a\) 的限制,那么我们要求的就是

\[\sum^{n}_{i=1} \sum^{m}_{j=1} \sigma(\gcd\left(i, j\right)) \]

可以转化为

\[\sum^{n}_{d=1} \sum^{n}_{i=1} \sum^{m}_{j=1} \sigma(d) \left[ \gcd\left(i, j\right)=d \right] \]

即

\[\sum^{n}_{d=1} \sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{i=1} \sum^{\left\lfloor \frac{m}{d}\right\rfloor}_{j=1} \sigma(d) \left[ \gcd\left(i, j\right)=1 \right] \]

利用常规的变换技巧

\[\sum^{n}_{d=1} \sigma(d) \sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{i=1} \sum^{\left\lfloor \frac{m}{d}\right\rfloor}_{j=1} \sum_{x \mid \gcd(i,j)}\mu(x) \]

换成整除分块形式

\[\sum^{n}_{d=1} \sigma(d) \sum^{\left\lfloor \frac{n}{d}\right\rfloor}_{x=1} \mu(x) {\left\lfloor \frac{n}{dx}\right\rfloor} {\left\lfloor \frac{m}{dx}\right\rfloor} \]

和上一道题类似,我们令 \(T=dx\),然后更改枚举顺序,有

\[\sum^{n}_{T=1} {\left\lfloor \frac{n}{T}\right\rfloor} {\left\lfloor \frac{m}{T}\right\rfloor} \sum_{d \mid T} \sigma(d) \mu(\frac{T}{d}) \]

这道题没有 \(a\) 的限制的部分就做完了,但还要考虑我们 \(10^9\) 的 \(a\)。

我们设 \(g(T)=\sum_{d \mid T}\sigma(d)\mu(\frac{T}{d})\),我们发现,只有在 \(\sigma(d) \leq a\) 时,才会对 \(g(T)\) 产生贡献。

那我们可以离线,把询问按照 \(a\) 从小到大排序,在从小到大逐个枚举询问的时候,\(a\) 不断变大,使得一些 \(\sigma(d)\) 对 \(g(T)\) 产生贡献。

我们枚举倍数来找到所有的 \(T\),需要动态修改 \(g(T)\) 且区间查询,考虑树状数组。

Code
#include <bits/stdc++.h>

using namespace std;

const int M = 1e5;
const int N = 100500;
const long long Mod = 1 << 31;

int T;

bool p[N];
int prime[N],cnt;
int num[N];
int sigma[N],mu[N];

unsigned long long ans[N];

struct Query{
    int id;
    int n,m,a;

    bool operator < (const Query &t) const {
        return a < t.a;
    }
}q[N];

struct Num{
    int x;

    bool operator < (const Num &t) const {
        return sigma[x] < sigma[t.x];
    }
}a[N];

void Seive() {
    sigma[1] = mu[1] = 1;
    for(int i = 2;i <= M; i++) {
        if(!p[i]) {
            cnt ++;
            prime[cnt] = i;
            mu[i] = -1;
            sigma[i] = num[i] = i + 1;
        }

        for(int j = 1;j <= cnt && i * prime[j] <= M; j++) {
            p[i * prime[j]] = 1;

            if(i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                num[i * prime[j]] = num[i] * prime[j] + 1;
                sigma[i * prime[j]] = sigma[i] / num[i] * num[i * prime[j]];

                break;
            }

            mu[i * prime[j]] = - mu[i];
            num[i * prime[j]] = prime[j] + 1;
            sigma[i * prime[j]] = sigma[i] * (prime[j] + 1);
        }
    }

    for(int i = 1;i <= M; i++)
        a[i].x = i;

    sort(a + 1,a + M + 1);
    return ;
}

int bit[N];

class BIT{
private:
    int lowbit(int x) {
        return x & (-x);
    }

public:
    void Add(int x,int y) {
        for(int i = x;i <= M; i += lowbit(i))
            bit[i] = (1ll * bit[i] + y) % Mod;
        
        return ;
    }

    int Query(int x) {
        int ans = 0;
        for(int i = x;i >= 1; i -= lowbit(i))
            ans = (1ll * ans + bit[i]) % Mod;
        
        return ans;
    }
}tree;

void Insert(int x) {
    for(int k = 1;x * k <= M; k++)
        tree.Add(x * k,1ll * mu[k] * sigma[x] % Mod);
        // 预处理
}

int query(int n,int m) {
    long long res = 0,last = 0,cur;

    for(int l = 1,r = 0;l <= n;l = r + 1,last = cur) {
        
        r = min(n / (n / l),m / (m / l));
        cur = tree.Query(r);
        res = (res + 1ll * (1ll * cur - last) * (n / l) % Mod * (m / l) % Mod) % Mod;
    }
    return res % Mod;
}

signed main() {
    ios::sync_with_stdio(false);
    cin >> T;

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

        if(q[i].n > q[i].m)
            swap(q[i].n,q[i].m);
    }

    sort(q + 1,q + T + 1);

    int j = 1;

    for(int i = 1;i <= T; i++) {
        for(    ;j <= M && sigma[a[j].x] <= q[i].a; j++) {
            Insert(a[j].x);
        }
        ans[q[i].id] = query(q[i].n,q[i].m);
    }

    for(int i = 1;i <= T; i++) 
        cout << (ans[i] % Mod + Mod) % Mod << "\n";
    return 0;
}

未完待续 \(\dots\)