首先有个东西是这题解题的关键：
[{
m{d}}(ij)=sum_{x|i}sum_{y|i}[{
m{gcd}}(x,y)=1]
]
有位dalao的题解证明了这个式子，可以去看看。
然后就可以开始推式子了：
[sum_{i=1}^nsum_{j=1}^md(ij) =sum_{i=1}^nsum_{j=1}^msum_{x|i}sum_{y|i}[{
m{gcd}}(x,y)=1]=sum_{x=1}^nsum_{y=1}^m[{
m{gcd}}(x,y)=1]lfloorfrac{n}{x}
floorlfloorfrac{m}{y}
floor
]
接下来用莫比乌斯函数替换掉 ([{ m{gcd}}(x,y)=1]) ：
[sum_{d|x,d|y}mu(d)=[{
m{gcd}}(x,y)=1]
]
代入：
[egin{align}
sum_{x=1}^nsum_{y=1}^m[{
m{gcd}}(x,y)=1]lfloorfrac{n}{x}
floorlfloorfrac{m}{y}
floor & =sum_{x=1}^nsum_{y=1}^m ig(sum_{d|x,d|y}mu(d)ig)lfloorfrac{n}{x}
floorlfloorfrac{m}{y}
floor \
end{align}
]
将这个式子由枚举 (x,y) 的形式转变为枚举 (d) ，即用 (xd,yd) 代替 (x,y) ：
[= sum_{d=1}^{{
m{min}}(n,m)}mu(d)sum_{x=1}^{lfloorfrac{n}{d}
floor}sum_{y=1}^{lfloorfrac{m}{d}
floor}lfloorfrac{n}{xd}
floorlfloorfrac{m}{yd}
floor
= sum_{d=1}^{{
m{min}}(n,m)}mu(d)sum_{x=1}^{lfloorfrac{n}{d}
floor}sum_{y=1}^{lfloorfrac{m}{d}
floor}lfloorfrac{lfloorfrac{n}{d}
floor}{x}
floorlfloorfrac{lfloorfrac{m}{d}
floor}{y}
floor
]
这里枚举的 (x,y) 是上式中的 (frac{x}{d},frac{y}{d}) 。
设 (f(n)=sum_{i=1}^{n}lfloorfrac{n}{i} floor) ，将其代入上式：
[sum_{i=1}^nsum_{j=1}^md(ij)=sum_{d=1}^{{
m{min}}(n,m)}mu(d)f(lfloorfrac{n}{d}
floor)f(lfloorfrac{m}{d}
floor)
]
到这里，很多题解就直接采用整除分块来计算 (f) 了，复杂度 (O(nsqrt{n})) 。但是其实还有一种更快的做法能做到 (O(n)) 预处理。
不难发现 (f(n)) 其实是计算了 (1 sim n) 每个数小于等于 (n) 的倍数个数之和，也就是说一个数对 (f(n)) 做出的贡献就是它的约数个数。那么有：
[f(n)=sum_{i=1}^{n}lfloorfrac{n}{i}
floor=sum_{i=1}^n{
m{d}}(i)
]
由于 ({ m{d}}) 是积性函数，可以使用线性筛做到 (O(n)) 预处理。
关于线性筛 ({ m{d}}) ，可以看这个。
({ m{Code}}:)
#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <cmath>
#define maxn 50005
#define Rint register int
#define INF 0x3f3f3f3f
using namespace std;
typedef long long lxl;
template <typename T>
inline T read()
{
T x=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}
while(ch>='0'&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return x*f;
}
int prime[maxn],cnt;
bool flag[maxn];
lxl d[maxn],f[maxn],mu[maxn],sum[maxn],num[maxn];
inline void sieve()
{
d[1]=1;
mu[1]=1;
for(int i=2;i<maxn;++i)
{
if(!flag[i]) prime[++cnt]=i,d[i]=2,num[i]=1,mu[i]=-1;
for(int j=1;j<=cnt&&i*prime[j]<maxn;++j)
{
flag[i*prime[j]]=true;
if(i%prime[j])
{
num[i*prime[j]]=1;
d[i*prime[j]]=d[i]*2;
}
else
{
num[i*prime[j]]=num[i]+1;
d[i*prime[j]]=d[i]/num[i*prime[j]]*(num[i*prime[j]]+1);
break;
}
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<maxn;++i)
f[i]=f[i-1]+d[i],sum[i]=sum[i-1]+mu[i];
}
inline lxl calcu(int n,int m)
{
lxl res=0;
if(n>m) swap(n,m);
for(int l=1,r=0;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
res+=f[n/l]*f[m/l]*(sum[r]-sum[l-1]);
}
return res;
}
int main()
{
// freopen("P3327.in","r",stdin);
sieve();
int T=read<int >();
while(T--)
{
int n=read<int >(),m=read<int >();
printf("%lld
",calcu(n,m));
}
return 0;
}