[Luogu P3327] [SDOI2015]约数个数和 (莫比乌斯反演)

题面：

传送门：洛咕

---

Solution

首先，我们需要一个结论：
\(\large d(i,j)=\sum_{x|i}\sum_{y|j}[gcd(x,y)=1]\)

证明
理性证明请看这篇博客的例五
本蒟蒻提供一个感性证明的方法：如果\(x*y\)是\(i*j\)的因数，我们必须有\(x|i,y|j\)，而后面那个\(gcd(x,y)\)是用来去重的

有了这个柿子之后，我们之后的推导就比较套路了：
为了方便讨论，之后的柿子均默认\(m>n\)
为了方便书写，之后的除法默认向下取整
原式：
\(\large \sum_{i=1}^n\sum_{j=1}^md(i*j)\)

把我们上面的结论代进去
\(\large \sum_{i=1}^n\sum_{j=1}^m\sum_{x|i}\sum_{y|j}[gcd(x,y)=1]\)

根据套路，这里的\(x|i\)与\(y|j\)应该写成枚举的形式：
\(\large \sum_{i=1}^n\sum_{j=1}^m\sum_{x=1}^{n}[x|i]\sum_{y=1}^{m}[y|j][gcd(x,y)=1]\)

这里显然可以把\(x,y\)的和式写到最前面去：
\(\large \sum_{x=1}^{n}[x|i]\sum_{y=1}^{m}[y|j]\sum_{i=1}^n\sum_{j=1}^m[gcd(x,y)=1]\)
然后就可以去掉\(x,y\)后面的那两个判断式啦
\(\large \sum_{x=1}^{n}\sum_{y=1}^{m}\sum_{i=1}^{n/x}\sum_{j=1}^{m/y}[gcd(x,y)=1]\)

然后我们再去套路后面那个\(gcd\)，根据莫比乌斯函数的性质\([x=1]=\sum_{d|x}\mu(d)\)，我们就把\(gcd\)带入得
\(\large \sum_{x=1}^{n}\sum_{y=1}^{m}\sum_{i=1}^{n/x}\sum_{j=1}^{m/y}\sum_{d|gcd(x,y)}\mu(d)\)
然后去枚举\(d\)
\(\large \sum_{x=1}^{n}\sum_{y=1}^{m}\sum_{i=1}^{n/x}\sum_{j=1}^{m/y}\sum_{d=1}^{n}\mu(d)[d|gcd(x,y)]\)
套路地把\(\mu\)的和式丢到最前面，化简一下就有：
\(\large \sum_{d=1}^{n}\mu(d)\sum_{x=1}^{n/d}\sum_{y=1}^{m/d}\sum_{i=1}^{n/(x*d)}\sum_{j=1}^{m/(y*d)}1\)
然后有
\(\large \sum_{d=1}^{n}\mu(d)\sum_{x=1}^{n/d}\sum_{y=1}^{m/d}\frac{n}{x*d}\frac{m}{y*d}\)

移项整理一下：
\(\large \sum_{d=1}^{n}\mu(d)\sum_{x=1}^{n/d}\frac{n}{x*d}\sum_{y=1}^{m/d}\frac{m}{y*d}\)

好了，到这里我就不会推了

之后的内容感谢@Maxwei_wzj的教学
事实上，这个柿子我们已经可以算了。
\(\large \sum_{d=1}^{n}\mu(d)\sum_{x=1}^{n/d}\frac{n}{x*d}\sum_{y=1}^{m/d}\frac{m}{y*d}\)
首先，我们有一个结论(我不会证)
\(\large x/(y*z)=(x/y)/z\) (这里的除法向下取整)

我们的柿子就可以变为
\(\large \sum_{d=1}^{n}\mu(d)\sum_{x=1}^{n/d}\frac{\frac{n}{d}}{x}\sum_{y=1}^{m/d}\frac{\frac{m}{d}}{y}\)

然后我们再设一个方程:
\(f(x)=\sum_{i=1}^x\frac{x}{i}\)
我们柿子就可以写为：
\(\large \sum_{d=1}^{n}\mu(d)f(\frac{n}{d})f(\frac{m}{d})\)

诶？我们好像又能整除分块了。
没错，我们只需要先用\(O(n*\sqrt n)\)的整除分块预处理出\(f\)
然后再每次\(O(\sqrt n)\)整除分块算出那个柿子就好。

时间复杂度\(O(n*\sqrt n)\)
完结撒花✿✿ヽ(°▽°)ノ✿

---

Code

//Luogu  P3327 [SDOI2015]约数个数和
//Jan,22ed,2019
//莫比乌斯反演
#include<iostream>
#include<cstdio>
using namespace std;
long long read()
{
	long long x=0,f=1; char c=getchar();
	while(!isdigit(c)){if(c=='-') f=-1;c=getchar();}
	while(isdigit(c)){x=x*10+c-'0';c=getchar();}
	return x*f;
}
const int N=50000+1000;
const int M=50000;
int prime[N],p_cnt,mu[N];
bool noPrime[N];
void GetPrime(int n)
{
	noPrime[1]=true,mu[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(noPrime[i]==false)
			prime[++p_cnt]=i,mu[i]=-1;
		for(int j=1;j<=p_cnt and i*prime[j]<=n;j++)
		{
			noPrime[i*prime[j]]=true;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=mu[i]*mu[prime[j]];
		}
	}
}
long long f[N],pre_mu[N];
int main()
{
	GetPrime(M);
	for(int i=1;i<=M;i++)
		for(int l=1,r;l<=i;l=r+1)
		{
			r=i/(i/l);
			f[i]+=(r-l+1)*(i/l);
		}
	for(int i=1;i<=M;i++)
		pre_mu[i]=pre_mu[i-1]+mu[i];
	
	int T=read();
	for(;T>0;T--)
	{
		long long n=read(),m=read();
		if(n>m) swap(n,m);
		
		long long ans=0;
		for(int l=1,r;l<=n;l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			ans+=(pre_mu[r]-pre_mu[l-1])*f[n/l]*f[m/l];
		}
		
		printf("%lld\n",ans);
	}
	return 0;
}