浅谈算法——莫比乌斯反演

前言

莫比乌斯反演(又称懵逼钨丝繁衍)，那种让人看了就懵逼的东西(其实是我太菜了)
莫比乌斯反演在知道之后对解题十分有帮助，\(O(n)\)的柿子分分钟化成\(O(\sqrt n)\)
那么，什么是莫比乌斯反演呢？

莫比乌斯反演

1.莫比乌斯反演
如果说，有\(f(n)\)和\(g(n)\)是定义在正整数集合上的两个函数，并且满足

\[f(n)=\sum\limits_{d|n} g(d) \]

那么就有

\[g(n)=\sum\limits_{d|n} f(\dfrac{n}{d})\times \mu(d) \]

这个变形也就是莫比乌斯反演

2.莫比乌斯函数(性质/证明)
对于任意整数n，满足

\[\sum\limits_{d|n} \mu(d)=\begin{cases}1&,n=1\\0&,n>1\end{cases} \]

证明：

• n=1时显然
• n>1时，设\(n=P_1^{x_1}\times P_2^{x_2}\times ...\times P_k^{x_k}\)，其中\(P_i(1\leqslant i\leqslant k)\)为n的互质因子，我们设\(d=P_1^{y_1}\times P_2^{y_2}\times ...\times P_k^{y_k}\)，其中\(0\leqslant y_i\leqslant x_i(1\leqslant i\leqslant k)\)，当存在一个\(y_i\geqslant2\)的话，那么\(\mu(d)=0\)，因此我们只要考虑\(y_i=0,1\)的情况。假设d中含有r个n的互质因子，且这些质因子的指数为1，那么\(\mu(d)=(-1)^r\)，这样的d一共有\(\binom{k}{r}\)个。
  那么\(\sum\limits_{d|n} \mu(d)=\sum\limits_{r=0}^k \binom{k}{r}\times (-1)^r=(1-1)^k=0\)

---

莫比乌斯函数\(\mu(n)\)是个积性函数，证明略。
这说明\(\mu(n)=\prod\limits_{i=1}^k \mu(P_i^{x_i})\)

---

由于莫比乌斯函数的这些性质，使得它可以用欧拉筛在线性时间内求出

3.莫比乌斯反演(性质/证明)
首先我们证明一下第一个柿子的正确性

\[g(n)=\sum\limits_{d|n}f(\dfrac{n}{d})\times \mu(d)=\sum\limits_{d|n}\mu(d) \sum\limits_{i|\frac{n}{d}}g(i)\Longrightarrow\sum\limits_{i|n}g(i)\sum\limits_{i\times d|n}\mu(d)=\sum\limits_{i|n}g(i)\sum\limits_{d|\frac{n}{i}}\mu(d) \]

箭头转换的时候，将\(g(i)\)提前枚举了
我们现在考虑一下\(\sum\limits_{d|\frac{n}{i}}\mu(d)\)的取值

• 当i=n时，\(\sum\limits_{d|\frac{n}{i}}\mu(d)=1,g(i)\sum\limits_{d|\frac{n}{i}}\mu(d)=g(n)\)
• 当i为小于n的约数时，\(\sum\limits_{d|\frac{n}{i}}\mu(d)=0,g(i)\sum\limits_{d|\frac{n}{i}}\mu(d)=0\)

所以\(g(n)=\sum\limits_{d|n}f(\dfrac{n}{d})\times \mu(d)\)得证

---

莫比乌斯函数有个性质：
当\(f(n)/g(n)\)中的任意一个为积性函数后，另一个也同样为积性函数
这个性质保证了其在之后能够进行线筛。
证明略(证明太长，我懒得写)

4.莫比乌斯函数的应用

• \(g(i)\)很难直接求，但\(\sum\limits_{d|i}g(d)\)很好求，我们就可以用莫比乌斯反演来求\(g(i)\):$$f(i)=\sum\limits_{d|i}g(d)\Rightarrow g(i)=\sum\limits_{d|i}f(\dfrac{i}{d})\times \mu(d)$$
• \(g(i)\)很难直接求，但\(\sum\limits_{d=1}^{\frac{n}{i}}g(d)\)很好求，我们就可以用莫比乌斯反演来求\(g(i)\):$$f(i)=\sum\limits_{d=1}^{\frac{n}{i}}g(d\times i)\Rightarrow g(i)=\sum\limits_{d=1}^{\frac{n}{i}}f(d\times i)\mu(d)$$

例题

[POI2007]Zap
Description
FGD正在破解一段密码，他需要回答很多类似的问题：对于给定的整数a,b和d，有多少正整数对x,y，满足x<=a，y<=b，并且gcd(x,y)=d。作为FGD的同学，FGD希望得到你的帮助。

Input
第一行包含一个正整数n，表示一共有n组询问。（1<=n<= 50000）接下来n行，每行表示一个询问，每行三个正整数，分别为a,b,d。（1<=d<=a,b<=50000）

Output
对于每组询问，输出到输出文件zap.out一个正整数，表示满足条件的整数对数。

Sample Input
2
4 5 2
6 4 3

Sample Output
3
2

HINT
对于第一组询问，满足条件的整数对有(2,2),（2,4）,（4,2）。对于第二组询问，满足条件的整数对有（6,3）,（3,3）。

Attention:为了方便书写，我们把x,y改成i,j，把多组数据的n改为T，把a,b改成n,m，并且我们让n<m

首先我们列出本题要求的柿子

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

d看起来太讨厌，我们把它直接除掉

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

然后我们发现最后面那一个长得和莫比乌斯函数很像，稍微改一下得到

\[\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum\limits_{x|i,x|j}\mu(x) \]

把枚举x的挪到前面

\[\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)\lfloor\dfrac{n}{dx}\rfloor\lfloor\dfrac{m}{dx}\rfloor \]

这个柿子终于变得小清新了，但是如果直接枚举x的话还是摆脱不了TLE的命运。。。
注意到后面的\(\lfloor\dfrac{n}{dx}\rfloor\lfloor\dfrac{m}{dx}\rfloor​\)在很多情况下值都是一样的，那么我们进行数论分块即可
总时间复杂度降为\(O(T\sqrt n)\)

/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline int read(){
	int x=0,f=1;char ch=getchar();
	for (;ch<'0'||ch>'9';ch=getchar())	if (ch=='-')    f=-1;
	for (;ch>='0'&&ch<='9';ch=getchar())	x=(x<<1)+(x<<3)+ch-'0';
	return x*f;
}
inline void print(int x){
	if (x>=10)     print(x/10);
	putchar(x%10+'0');
}
const int N=5e4;
int prime[N+10],miu[N+10],sum[N+10];
bool inprime[N+10];
void prepare(){
	miu[1]=1;
	int tot=0;
	for (int i=2;i<=N;i++){
		if (!inprime[i])	prime[++tot]=i,miu[i]=-1;
		for (int j=1;j<=tot&&i*prime[j]<=N;j++){
			inprime[i*prime[j]]=1;
			if (i%prime[j]==0){miu[i*prime[j]]=0;break;}
			miu[i*prime[j]]=-miu[i];
		}
	}
	for (int i=1;i<=N;i++)	sum[i]=sum[i-1]+miu[i];
}
int main(){
	prepare();
	for (int Data=read();Data;Data--){
		int A=read(),B=read(),D=read();
		ll Ans=0;
		A/=D,B/=D;
		int x=min(A,B),pos=0;
		for (int d=1;d<=x;d=pos+1){
			pos=min(A/(A/d),B/(B/d));
			Ans+=1ll*(sum[pos]-sum[d-1])*(A/d)*(B/d);
		}
		printf("%lld\n",Ans);
	}
	return 0;
}