P3327 题解（莫反）
简要题意：
设 \(d(x)\) 为 \(x\) 的约数个数，给定 \(n,m\)，求：
\[\sum_{i=1}^n\sum_{j=1}^md(ij)
\]
多组测试数据
首先，我们可以证明：
\[d(ij)=\sum\limits_{x|i}\sum\limits_{y|j}[gcd(x,y)=1]
\]
考虑将数一一映射：若\(d(ij)\)中有约数\(p^c\)，\(i\)中\(p\)的最高次幂为\(p^a\)，\(j\)中\(p\)的最高次幂为\(p^b\)，必然有\(c\leqslant a+b\)，则我们定义规则：
若\(c<a\)，则我们在\(i\)中取\(p^c\)，在\(j\)中取\(1\)，即令\(x*p^c,y*1\)
若\(c>a\)，则我们在\(j\)中取\(p^{(c-a)}\)，在\(i\)中取\(1\)，即令\(x*1,y*p^{c-a}\)
我们发现，对于\(ij\)中的每个约数都可以映射出来，且满足\(gcd(x,y)=1\)，反之，满足\(gcd(x,y)=1\)的每个数，都对应的\(ij\)的一个约数
证毕
则原式可表示成：
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{x|i}\sum\limits_{y|j}[gcd(x,y)=1]
\]
将\(x,y\)提到前面来：
\[\sum\limits_{x=1}^n\sum\limits_{y=1}^m[gcd(x,y)=1]\left\lfloor\frac{n}{x}\right\rfloor\left\lfloor\frac{m}{y}\right\rfloor
\]
令：
\[f(d)=\sum\limits_{x=1}^{n}\sum\limits_{y=1}^{m}[gcd(x,y)=d]\left\lfloor\frac{n}{x}\right\rfloor\left\lfloor\frac{m}{y}\right\rfloor
\]
\[g(d)=\sum\limits_{d|i}f(i)=\sum\limits_{x=1}^{n/d}\sum\limits_{y=1}^{n/d}\left\lfloor\frac{n}{dx}\right\rfloor\left\lfloor\frac{m}{dy}\right\rfloor
\]
则原式莫反一下为为：
\[f(1)=\sum\limits_{1|i}\mu(i/1)g(i)=\sum\limits_{i=1}^n\mu(i)g(i)=\sum\limits_{i=1}^n\mu(i)\sum\limits_{x=1}^{n/i}\sum\limits_{y=1}^{n/i}\left\lfloor\frac{n}{ix}\right\rfloor\left\lfloor\frac{m}{iy}\right\rfloor
\]
预处理一下，令：
\[s(d)=\sum\limits_{i=1}^d\left\lfloor\frac{d}{i}\right\rfloor
\]
则有：
\[g(d)=s(\left\lfloor\frac{n}{d}\right\rfloor)*s(\left\lfloor\frac{m}{d}\right\rfloor)
\]
\[f(1)=\sum\limits_{i=1}^n\mu(i)s(\left\lfloor\frac{n}{d}\right\rfloor)s(\left\lfloor\frac{m}{d}\right\rfloor)
\]
分块计算\(s(\left\lfloor\frac{n}{d}\right\rfloor)s(\left\lfloor\frac{m}{d}\right\rfloor)\)即可
上代码：
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const ll N=5e4+50,mx=5e4;
ll T,n,m;
ll zs[N],mu[N],s[N],sm[N],idx,ans;
bool prime[N];
void init()
{
mu[1]=1;
for(ll i=2;i<=mx;i++)
{
if(!prime[i]) zs[++idx]=i,mu[i]=-1;
for(ll j=1;j<=idx;j++)
{
ll x=zs[j]*i;
if(x>mx) break;
prime[x]=true;
if(i%zs[j]==0) break;
mu[x]=-mu[i];
}
}
for(ll i=1;i<=mx;i++)
{
sm[i]=sm[i-1]+mu[i];
for(ll l=1,r;l<=i;l=r+1)
{
r=i/(i/l);
s[i]=s[i]+(i/l)*(r-l+1);
}
}
}
int main()
{
init();
scanf("%lld",&T);
while(T--)
{
scanf("%lld %lld",&n,&m);
if(n>m) swap(n,m);
ans=0;
for(ll l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans=ans+(sm[r]-sm[l-1])*s[n/l]*s[m/l];
}
printf("%lld\n",ans);
}
return 0;
}