莫比乌斯函数及莫比乌斯反演

莫比乌斯函数

\[\mu(x)=\begin{cases}1\ \ &(x=1)\\(-1)^k\ \ &(\text{$x$没有平方数因子，且$x$的质因子个数为$k$})\\0 &(\text{$x$有平方数因子})\end{cases} \]

很好理解 \(\mu\) 是积性函数。但是这个定义略显奇怪。这是因为它满足如下式子：

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

也就是说有

\[\mu * 1=\varepsilon \]

事实上莫比乌斯函数是为了满足上面的那个式子才被定义出来的。

类似的有

\[\sum_{d|\gcd(i,j)}\mu(d)=[\gcd(i,j)=1] \]

事实上只是将 \(n\) 换成了 \(\gcd(i,j)\)，但是在很多式子中会很有用。

莫比乌斯反演

\[g(n)=\sum_{i|n}f(i)\Leftrightarrow f(n)=\sum_{i|n}\mu(\frac{n}{i})g(i) \]

但是事实上这个式子并不是很常用，限制有点多。
更常用的是上面的式子和另一个，把它们合起来写：

\[\begin{cases} \sum_{d|n}\mu(d)=[n=1] \\ n=\sum_{d|n} \varphi(d) \end{cases} \]

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

求

\[\sum_{i=1}^n\sum_{j=1}^m\operatorname{lcm}(i, j) \]

显然先将 \(\operatorname{lcm}\) 拆开，拆成 \(\gcd\) 后套路性的枚举 \(\gcd\)，类似交换和式处理使得出现整除分块的形式：

\[=\sum_{d=1}^nd\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[\gcd(i,j)=1]ij \]

我们把 \(\sum_{i=1}^{x}\sum_{j=1}^{y}[\gcd(i,j)=1]ij\) 设为 \(solve1(x,y)\)，那上面的式子就是

\[\sum_{d=1}^nd \ solve1(\left \lfloor n/d \right \rfloor,\left \lfloor m/d \right \rfloor ) \]

显然可以整除分块。然后再来考虑 \(solve1\) 怎么处理。

用上面的式子把 \([\gcd(i,j)=1]\) 拆成 \(\sum_{d|\gcd(i,j)}\mu(d)\)，仍然是交换和式变成可整除分块的形式：

\[\begin{aligned} solve1(n,m)&=\sum_{i=1}^n\sum_{j=1}^mij\sum_{k|\gcd(i,j)}\mu(k)\\ &=\sum_{k=1}^n\mu(k)k^2\sum_{i=1}^{n/k}i\sum_{j=1}^{m/k}j \end{aligned} \]

这时最后两个求和符号显然都可以 \(O(1)\) 计算。记 \(g(n,m)=\sum_{i=1}^n i \sum_{j=1}^n j=\frac{i(i+1)}{2}\frac{j(j+1)}{2}\)，

\[=\sum_{k=1}^n\mu(k)k^2g(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor) \]

这样 \(solve1\) 也可以整除分块。两个整除分块的嵌套听大佬说时间复杂度是 \(O(n)\) 的。（不太会证）

code

点击查看代码

#include<bits/stdc++.h>
using namespace std;
#define int long long 
const int N=1e7+7;
const int p=20101009;
int n,m,mu[N],zhi[N],sum[N],ni;
bool sign[N];
int ksm(int x,int k){int res=1;while(k){if(k&1)res=res*x%p;x=x*x%p,k>>=1;}return res;}
void init()
{
	mu[1]=1;int cnt=0,up=min(n,m);sign[0]=sign[1]=1;ni=ksm(2,p-2);
	for(int i=2;i<=up;i++){
		if(!sign[i]) zhi[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*zhi[j]<=up;j++){
			sign[i*zhi[j]]=1;
			if(i%zhi[j]==0) {mu[i*zhi[j]]=0;break;}
			mu[i*zhi[j]]=-mu[i];
		}
	}
	for(int i=1;i<=up;i++) sum[i]=(sum[i-1]+i*i%p*mu[i]%p+p)%p;
}
int g(int x,int y){return x*(x+1)%p*y%p*(y+1)%p*ni%p*ni%p;}
int solve2(int x,int y)
{
	int res=0;
	for(int l=1,r;l<=min(x,y);l=r+1)
	{
		r=min(x/(x/l),y/(y/l));
		res=(res+(sum[r]-sum[l-1])%p*g(x/l,y/l)%p+p)%p;
	}
	return res;
}
int solve1()
{
	int res=0;
	for(int l=1,r;l<=min(n,m);l=r+1)
	{
		r=min(n/(n/l),m/(m/l));
		res=(res+(l+r)*(r-l+1)%p*ni%p*solve2(n/l,m/l)%p+p)%p;
	}
	return res;
}
signed main()
{
	cin>>n>>m;init();
	cout<<solve1()<<'\n';return 0;
}