莫比乌斯反演学习笔记+[POI2007]Zap（洛谷P3455，BZOJ1101）

先看一道例题：
[POI2007]Zap

BZOJ

洛谷

题目大意：$T$ 组数据，求 $sum^n_{i=1}sum^m_{j=1}[gcd(i,j)=k]$

$1leq Tleq 50000,1leq kleq n,mleq 50000$

暴力做法 $O(Tnmlogmax(n,m))$ 不用说了，那有没有什么更好的做法呢？

---

我们定义一种函数叫莫比乌斯函数 $mu$，它的定义是：

当 $n=1$ 时，$mu(n)=1$

当 $n$ 可以分解成 $p_1p_2...p_k$，其中 $p_i$ 均为质数且互不相同时，$mu(n)=(-1)^k$

否则 $mu(n)=0$

莫比乌斯函数有一个性质是这样的：

$sum_{d|n}mu(d)=[n=1]$

这个性质暂时用不到，以后到狄利克雷卷积和杜教筛时才有用，所以只做了解。

---

开始正文：莫比乌斯反演。

情况一：

若函数 $F$ 和 $f$ 满足 $F(n)=sum_{d|n}f(d)$

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

情况二：（用的较多，一定要死记硬背）

若函数 $F$ 和 $f$ 满足 $F(n)=sum_{n|d}f(d)$

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

让我们来感性理解一下。设 $F$ 和 $f$ 满足：

$F(1)=f(1)$

$F(2)=f(1)+f(2)$

$F(3)=f(1)+f(3)$

$F(4)=f(1)+f(2)+f(4)$

$F(5)=f(1)+f(5)$

$F(6)=f(1)+f(2)+f(3)+f(6)$

$F(7)=f(1)+f(7)$

$F(8)=f(1)+f(2)+f(4)+f(8)$

$F(9)=f(1)+f(3)+f(9)$

$F(10)=f(1)+f(2)+f(5)+f(10)$

$...$

那么就有：

$f(1)=F(1)$

$f(2)=F(2)-F(1)$

$f(3)=F(3)-F(1)$

$f(4)=F(4)-F(2)$

$f(5)=F(5)-F(1)$

$f(6)=F(6)-F(3)-F(2)+F(1)$

$f(7)=F(7)-F(1)$

$f(8)=F(8)-F(4)$

$f(9)=F(9)-F(3)$

$f(10)=F(10)-F(5)-F(2)+F(1)$

我们单独把 $f(6)$ 提出来看（其他的类似）

$f(6)=F(6)-F(3)-F(2)+F(1)=mu(1)F(6)+mu(2)F(3)+mu(3)F(2)+mu(6)F(1)=sum_{d|6}mu(d)F(frac{6}{d})$

实际上就是一个容斥原理。

证明较难，貌似要用狄利克雷卷积，此处略去。（其实是因为我太蒻了）

另外 $mu$ 还是个积性函数，虽然现在也没什么用。

如何线性筛 $mu$？

我们发现 $mu(1)=1,mu(prime)=-1$。

线筛的原理是搜到重复的质因子时退出，正好符合 $mu$ 的第三条！

所以：

当 $prime[j]|i$ 时，$mu(i imes prime[j])=0$

否则，$mu(i imes prime[j])=-mu(i)$

程序如下：

void init(int n){
    memset(vis,0,sizeof(vis));
    memset(mu,0,sizeof(mu));
    memset(prime,0,sizeof(prime));
    len=0;
    vis[1]=true;
    mu[1]=1;    //预处理1
    for(int i=2;i<=n;i++){
        if(!vis[i]){
            mu[i]=-1;    //是质数，莫比乌斯函数=-1
            prime[++len]=i;
        }
        for(int j=1;j<=len && i*prime[j]<=n;j++){
            int k=i*prime[j];
            vis[k]=true;
            if(i%prime[j]==0) break;    //有重复质因子，莫比乌斯函数=0
            else mu[k]=-mu[i];    //多了一个不重复质因子，莫比乌斯函数区相反数
        }
    }
}

莫比乌斯反演大部分题目都含有 $gcd$，套路就看例题，大部分题目都一样的。

---

回到原题。 

$T$ 组数据，求 $sum^n_{i=1}sum^m_{j=1}[gcd(i,j)=k]$

开始讲这类题目的套路：

假设 $n<m$。

设两个函数：

$f(x)=sum^n_{i=1}sum^m_{j=1}[gcd(i,j)=x]$

$F(x)=sum^n_{i=1}sum^m_{j=1}[x|gcd(i,j)]$

题目要求即为 $f(k)$。

我们发现在 $F(x)$ 中有序对 $(i,j)$ 对答案作出 $1$ 的贡献当且仅当 $x|i$ 且 $x|j$。

这样的 $i$ 有 $lfloorfrac{n}{x} floor$ 个，$j$ 有 $lfloorfrac{m}{x} floor$ 个。

所以 $F(x)=lfloorfrac{n}{x} floorlfloorfrac{m}{x} floor$

根据定义，$F(x)=sum^n_{x|d}f(d)$

莫比乌斯反演一波：$f(x)=sum^n_{x|d}mu(frac{d}{x})F(d)$

题目要求就变成了：$f(k)=sum^n_{k|d}mu(frac{d}{k})F(d)$

我们发现当且仅当 $d$ 是 $k$ 的倍数时对答案有贡献，那我们可以改一下枚举的方式：

$f(k)=sum^{lfloorfrac{n}{k} floor}_{d=1}mu(d)F(dk)$

把 $F(dk)$ 替换：$f(k)=sum^{lfloorfrac{n}{k} floor}_{d=1}mu(d)lfloorfrac{n}{dk} floorlfloorfrac{m}{dk} floor$

$d$ 看着不爽：$f(k)=sum^{lfloorfrac{n}{k} floor}_{i=1}mu(i)lfloorfrac{n}{ik} floorlfloorfrac{m}{ik} floor$

此时这个式子已经可以做到 $O(n)$ 计算了，线性筛出 $mu$ 然后扫一遍就行了。

等等，$T$ 组数据，$O(Tn)$？

---

阅读以下内容以前请先学会前置技能整除分块

 我们发现这里有个很明显的整除分块的形式，那么我们可以考虑 $[l,r]$ 这段区间，其中 $lfloorfrac{n}{ik} floor=lfloorfrac{n}{jk} floor=x$ 且 $lfloorfrac{m}{ik} floor=lfloorfrac{m}{jk} floor=y:i,jin[l,r]$

$ sum^r_{i=l}mu(i)lfloorfrac{n}{ik} floorlfloorfrac{m}{ik} floor$

$=sum^r_{i=l}mu(i)xy$

$=xysum^r_{i=l}mu(i)$

那么我们只需要求出 $mu$ 的前缀和，然后整除分块套上去即可。

还可以加一个常数优化：

我们发现，在原式中，只要出现了 $n$ 和 $m$ 的地方都是 $lfloorfrac{n}{ik} floor$ 和 $lfloorfrac{m}{ik} floor$ 的形式。

考虑到 $lfloorfrac{n}{ik} floor=lfloorfrac{frac{n}{k}}{i} floor$，我们可以在开始整除分块之前就 $n$ 和 $m$ 除以 $k$ 然后再分块，可以少掉一个 $sqrt{k}$ 的常数。

---

代码如下：时间复杂度 $O(Tsqrt{n})$，空间复杂度 $O(n)$

#include<bits/stdc++.h>
using namespace std;
int t,n,m,k;
int prime[50050],mu[50050],pre[50050],len;
bool vis[50050];
void init(int x){
    vis[1]=true;
    mu[1]=1;
    for(int i=2;i<=x;i++){
        if(!vis[i]){
            mu[i]=-1;
            prime[++len]=i;
        }
        for(int j=1;j<=len && i*prime[j]<=x;j++){
            int k=i*prime[j];
            vis[k]=true;
            if(i%prime[j]==0) break;
            else mu[k]=-mu[i];
        }
    }
    for(int i=1;i<=x;i++) pre[i]=pre[i-1]+mu[i];
}
int main(){
    init(50000);
    scanf("%d",&t);
    while(t--){
        scanf("%d%d%d",&n,&m,&k);
        n/=k;m/=k;
        int ans=0;
        for(int l=1,r;l<=min(n,m);l=r+1){
            r=min(n/(n/l),m/(m/l));
            ans+=(n/l)*(m/l)*(pre[r]-pre[l-1]);
        }
        printf("%d
",ans);
    }
}

然后推荐几题：

洛谷P2522  BZOJ2301 [HAOI2011]Problem b （题解待填充）

洛谷P2257 YY的GCD （题解待填充）

洛谷P1447 BZOJ2005 [NOI2010]能量采集 （题解待填充）