我的做法基于以下两个公式：
[[n=1]=sum_{d|n}mu(d)
]
[sigma_0(i*j)=sum_{x|i}sum_{y|j}[gcd(x,y)=1]
]
其中(sigma_0(n))表示(n)的约数个数
第一个公式是莫比乌斯函数的基本性质，至于第二个公式的证明，可以考虑(i*j)中每一个质因子对 (sigma_0(i*j)) 的贡献，对于一个质因子 (p) ，若它在 (i) 中的次数为 (k_1) ，它在 (j) 中的次数为 (k_2) ，则在 (sigma_0(i*j)) 中(p)的贡献为((k_1+k_2+1))（约数个数计算公式），而在(gcd(x,y)=1)的情况下，要么(x)中(p)的次数为0，要么(y)中(p)的次数为0，一共有((k_1+k_2+1))种方案，与(i*j)中的贡献相同，所以等式左右两边相等。
然后就可以愉快的推式子啦！
[sum_{i=1}^{n} sum_{j=1}^{m}sigma_0(i*j)
]
[sum_{i=1}^{n} sum_{j=1}^{m}sum_{x|i}sum_{y|j}[gcd(x,y)=1]
]
[sum_{x=1}^{n}sum_{y=1}^m sum_{i=1}^{lfloor frac{n}{x}
floor } sum_{j=1}^{lfloor frac{m}{y}
floor} [gcd(x,y)=1]
]
[sum_{x=1}^{n}sum_{y=1}^m lfloor frac{n}{x}
floor lfloor frac{m}{y}
floor sum_{k|gcd(x,y)}mu(k)
]
[sum_{k=1}^{n}sum_{x=1}^{lfloor frac{n}{k}
floor}sum_{y=1}^{lfloor frac{m}{k}
floor} lfloor frac{n}{kx}
floor lfloor frac{m}{ky}
floormu(k)
]
[sum_{k=1}^{n}mu(k)(sum_{x=1}^{lfloor frac{n}{k}
floor}lfloorfrac{n}{kx}
floor)(sum_{y=1}^{lfloorfrac{m}{k}
floor}lfloorfrac{m}{ky}
floor)
]
然后我们设(S(n)=sum_{i=1}^{n}lfloorfrac{n}{i} floor)，显然(S(n))是可以(O(sqrt{n}))计算的
则上式可化为：
[sum_{k=1}^{n}mu(k)S(lfloorfrac{n}{k}
floor)S({lfloorfrac{m}{k}
floor})
]
先预处理(S(1)-S(maxn))，然后就可以(O(sqrt{n}))回答每组询问啦！
代码：
#include<bits/stdc++.h>
using namespace std;
#define N 50007
#define ll long long
const int lim=50000;
ll s[N];
int ui[N],pr[N],cnt;
bool zhi[N];
void Init()
{
int i,j;
ui[1]=1;
for(i=2;i<=lim;i++)
{
if(!zhi[i])
{
pr[++cnt]=i;
ui[i]=-1;
}
for(j=1;j<=cnt&&i*pr[j]<=lim;j++)
{
int p=pr[j],x=i*p;
zhi[x]=true;
if(i%p==0)
{
ui[x]=0;
break;
}
ui[x]=-ui[i];
}
}
for(i=1;i<=lim;i++)
ui[i]+=ui[i-1];
for(i=1;i<=lim;i++)
{
int l,r;
for(l=1;l<=i;l=r+1)
{
r=i/(i/l);
s[i]+=1ll*(r-l+1)*(i/l);
}
}
}
int main()
{
int n,m,t;
Init();
scanf("%d",&t);
while(t--)
{
scanf("%d%d",&n,&m);
int l1=1,r1,l2=1,r2,cur=1;
ll ans=0;
while(l1<=n&&l2<=m)
{
int l,r;
r1=n/(n/l1),r2=m/(m/l2);
l=cur;
if(r1<r2)r=r1,cur=l1=r1+1;
else r=r2,cur=l2=r2+1;
ans+=1ll*(ui[r]-ui[l-1])*s[n/l]*s[m/l];
}
printf("%lld
",ans);
}
return 0;
}