莫比乌斯反演小记

基本内容

莫比乌斯函数

\(\mu\) 定义为 \(1\) 的逆。

一些小性质：

1. \(\mu * 1=\epsilon\)
2. \(\mu * \text{id}=\varphi\)

反演内容

我的理解是：

\[[a=1]=\sum\limits_{d|a}\mu(d) \]

典型例题

例1 P2398 GCD SUM

求

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

来推下式子：

先枚举 \(\gcd\)：

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

再将 \(d\) 除掉：

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

然后莫比乌斯反演：

\[\sum\limits_{d=1}^n d \sum\limits_{i=1}^{\left\lfloor \frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor \frac{n}{d}\right\rfloor}\sum\limits_{k|\gcd(i,j)}\mu(k) \]

枚举 \(k\)：

\[\sum\limits_{k=1}^n\mu(k)\sum\limits_{d=1}^n d \sum\limits_{i=1}^{\left\lfloor \frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor \frac{n}{d}\right\rfloor}[k|\gcd(i,j)] \]

将 \(k\) 除掉：

\[\sum\limits_{k=1}^n\mu(k)\sum\limits_{d=1}^n d \sum\limits_{i=1}^{\left\lfloor \frac{n}{dk}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor \frac{n}{dk}\right\rfloor}*1 \]

令 \(T=dk\)，整理一下：

\[\sum\limits_{T=1}^n\left\lfloor \frac{n}{T}\right\rfloor^2\sum\limits_{k|T}\mu(k)(\frac{T}{k}) \]

由于 \(\mu*\text{id}=\varphi\)：

\[\sum\limits_{T=1}^n\left\lfloor \frac{n}{T}\right\rfloor^2\varphi(T) \]

欧拉函数可以直接线性筛预处理前缀和。中间的可以用整除分块 \(\Theta(\sqrt n)\)解决。

时间复杂度瓶颈在线性筛的 \(\Theta(n)\)。

点击查看代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN=2e5+5;

int n,cnt;
int prime[MAXN],vis[MAXN],phi[MAXN],sum[MAXN];

inline void init(int n)
{
    memset(vis,0,sizeof(vis));
    phi[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            prime[++cnt]=vis[i]=i;
            phi[i]=i-1;
        }
        for(int j=1;j<=cnt && i*prime[j]<=n;j++)
        {
            vis[i*prime[j]]=prime[j];
            if(i%prime[j]) phi[i*prime[j]]=phi[i]*(prime[j]-1);
            else phi[i*prime[j]]=phi[i]*prime[j];
        }
    }
    for(int i=1;i<=n;i++) sum[i]=sum[i-1]+phi[i];
    return;
}

signed main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(0),cout.tie(0);
    cin>>n;
    init(n);
    int ans=0;
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans+=(sum[r]-sum[l-1])*(n/l)*(n/l);
    }
    printf("%lld\n",ans);
    return 0;
}

例2 P1829 [国家集训队] Crash的数字表格 / JZPTAB

求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\text{lcm}(i,j) \]

这道题与上题差不多，因为 \(\text{lcm}(i,j)=\frac{i\times j}{\gcd(i,j)}\)。

推式子：

令 \(N=\min(n,m)\)，

\[\begin{aligned} \sum\limits_{i=1}^n\sum\limits_{j=1}^m\text{lcm}(i,j)&=\sum\limits_{i=1}^N\sum\limits_{j=1}^N\frac{ij}{\gcd(i,j)}\\ &=\sum\limits_{d=1}^N\sum\limits_{i=1}^N\sum\limits_{j=1}^N\frac{ij}{[\gcd(i,j)=d]}\\ &=\sum\limits_{d=1}^N\sum\limits_{i=1}^N\sum\limits_{j=1}^N ij[\gcd(i,j)=d]\\ &=\sum\limits_{d=1}^N d\sum\limits_{i=1}^{\left\lfloor \frac{N}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor \frac{N}{d}\right\rfloor} ij[\gcd(i,j)=1]\\ &=\sum\limits_{d=1}^N d\sum\limits_{i=1}^{\left\lfloor \frac{N}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor \frac{N}{d}\right\rfloor} ij\sum\limits_{k|\gcd(i,j)}\mu(k)\\ &=\sum\limits_{k=1}^N\mu(k)\sum\limits_{d=1}^N d\sum\limits_{i=1}^{\left\lfloor \frac{N}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor \frac{N}{d}\right\rfloor} ij[k|\gcd(i,j)]\\ &=\sum\limits_{k=1}^N\mu(k)k^2\sum\limits_{d=1}^N d\sum\limits_{i=1}^{\left\lfloor \frac{N}{dk}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor \frac{N}{dk}\right\rfloor} ij\\ \end{aligned} \]

然后到这里就可以做了，莫比乌斯函数预处理，然后两次整除分块，时间复杂度 \(\Theta(n)\)。

点击查看代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN=1e7+5;
const int MOD=20101009;

int n,m,cnt;
int prime[MAXN],vis[MAXN],mu[MAXN],sum[MAXN];

inline void init(int n)
{
    mu[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i]) {mu[i]=-1;prime[++cnt]=i;}
        for(int j=1;j<=cnt && i*prime[j]<=n;j++)
        {
            vis[i*prime[j]]=1;
            if(!(i%prime[j])) {mu[i*prime[j]]=0;break;}
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;i++) sum[i]=(sum[i-1]+i*i%MOD*(mu[i]+MOD)%MOD)%MOD;
    return;
}

inline int S(int r,int l=1) {return ((r+l)*(r-l+1)/2)%MOD;}
inline int Sum(int x,int y) {return ((x*(x+1)/2)%MOD)*((y*(y+1)/2)%MOD)%MOD;}

inline int solve(int x,int y)
{
    int ans=0;
    for(int l=1,r;l<=min(x,y);l=r+1)
    {
        int r1=x/(x/l),r2=y/(y/l);
        r=min(r1,r2);
        ans=(ans+(sum[r]-sum[l-1]+MOD)%MOD*Sum(x/r,y/r)%MOD)%MOD;
    }
    return ans;
}

signed main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(0),cout.tie(0);
    cin>>n>>m; init(min(n,m));
    int ans=0;
    for(int l=1,r;l<=min(n,m);l=r+1)
    {
        int r1=n/(n/l),r2=m/(m/l);
        r=min(r1,r2);
        ans=(ans+S(r,l)*solve(n/r,m/r)%MOD)%MOD;
    }
    printf("%lld\n",(ans+MOD)%MOD);
    return 0;
}

例3 P3327 [SDOI2015] 约数个数和

这道题有个前置知识：

\[d(ij)=\sum\limits_{x|i}\sum\limits_{y|j}[\gcd(x,y)=1] \]

然后就来推式子：

令 \(N=\min(n,m)\)，

\[\begin{aligned} \sum\limits_{i=1}^N\sum\limits_{j=1}^N\sum\limits_{x|i}\sum\limits_{y|j}[\gcd(i,j)=1]&=\sum\limits_{i=1}^N\sum\limits_{j=1}^N\left\lfloor\frac{n}{i}\right\rfloor \left\lfloor\frac{m}{j}\right\rfloor[\gcd(i,j)=1]\\ &=\sum\limits_{i=1}^N\sum\limits_{j=1}^N\left\lfloor\frac{n}{i}\right\rfloor \left\lfloor\frac{m}{j}\right\rfloor\sum\limits_{d|\gcd(i,j)}\mu(d )\\ &=\sum\limits_{d=1}^N\mu(d)\sum\limits_{i=1}^N\sum\limits_{j=1}^N\left\lfloor\frac{n}{i}\right\rfloor \left\lfloor\frac{m}{j}\right\rfloor[d|\gcd(i,j)]\\ &=\sum\limits_{d=1}^N\mu(d)\sum\limits_{i=1}^{\left\lfloor\frac{N}{i}\right\rfloor }\sum\limits_{j=1}^{\left\lfloor\frac{N}{i}\right\rfloor }\left\lfloor\frac{n}{id}\right\rfloor \left\lfloor\frac{m}{jd}\right\rfloor \end{aligned} \]

那么到这里就可以做了。

显然 \(s(x)=\sum\limits_{i=1}^x\left\lfloor\frac{x}{i}\right\rfloor\) 可以整除分块预处理，那么每组数据计算只需要用一次整除分块，整体时间复杂度 \(\Theta(T\sqrt n)\)。

点击查看代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN=5e4+5;

int T,n,m,cnt;
int prime[MAXN],vis[MAXN],mu[MAXN],summu[MAXN],sum[MAXN];

inline void init(int n)
{
    mu[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i]) {mu[i]=-1;prime[++cnt]=i;}
        for(int j=1;j<=cnt && i*prime[j]<=n;j++)
        {
            vis[i*prime[j]]=1;
            if(!(i%prime[j])) {mu[i*prime[j]];break;}
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;i++) summu[i]=summu[i-1]+mu[i];
    for(int i=1;i<=n;i++)
    {
        int ans=0;
        for(int l=1,r;l<=i;l=r+1)
            r=i/(i/l),ans+=(r-l+1)*(i/l);
        sum[i]=ans;
    }
    return;
}

signed main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(0),cout.tie(0);
    cin>>T; init(MAXN);
    while(T--)
    {
        cin>>n>>m; int ans=0;
        for(int l=1,r;l<=min(n,m);l=r+1)
        {
            int r1=n/(n/l),r2=m/(m/l);
            r=min(r1,r2);
            ans+=(summu[r]-summu[l-1])*sum[n/l]*sum[m/l];
        }
        printf("%lld\n",ans);
    }
    return 0;
}