首先有个东西是这题解题的关键：

\[{\rm{d}}(ij)=\sum_{x|i}\sum_{y|i}[{\rm{gcd}}(x,y)=1] \]

有位dalao的题解证明了这个式子，可以去看看。

然后就可以开始推式子了：

\[\sum_{i=1}^n\sum_{j=1}^md(ij) =\sum_{i=1}^n\sum_{j=1}^m\sum_{x|i}\sum_{y|i}[{\rm{gcd}}(x,y)=1]=\sum_{x=1}^n\sum_{y=1}^m[{\rm{gcd}}(x,y)=1]\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor \]

接下来用莫比乌斯函数替换掉 \([{\rm{gcd}}(x,y)=1]\) ：

\[\sum_{d|x,d|y}\mu(d)=[{\rm{gcd}}(x,y)=1] \]

代入：

\[\begin{align} \sum_{x=1}^n\sum_{y=1}^m[{\rm{gcd}}(x,y)=1]\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor & =\sum_{x=1}^n\sum_{y=1}^m \big(\sum_{d|x,d|y}\mu(d)\big)\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor \\ \end{align} \]

将这个式子由枚举 \(x,y\) 的形式转变为枚举 \(d\) ，即用 \(xd,yd\) 代替 \(x,y\) ：

\[= \sum_{d=1}^{{\rm{min}}(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{n}{xd}\rfloor\lfloor\frac{m}{yd}\rfloor = \sum_{d=1}^{{\rm{min}}(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{x}\rfloor\lfloor\frac{\lfloor\frac{m}{d}\rfloor}{y}\rfloor \]

这里枚举的 \(x,y\) 是上式中的 \(\frac{x}{d},\frac{y}{d}\) 。

设 \(f(n)=\sum_{i=1}^{n}\lfloor\frac{n}{i}\rfloor\) ，将其代入上式：

\[\sum_{i=1}^n\sum_{j=1}^md(ij)=\sum_{d=1}^{{\rm{min}}(n,m)}\mu(d)f(\lfloor\frac{n}{d}\rfloor)f(\lfloor\frac{m}{d}\rfloor) \]

到这里，很多题解就直接采用整除分块来计算 \(f\) 了，复杂度 \(O(n\sqrt{n})\) 。但是其实还有一种更快的做法能做到 \(O(n)\) 预处理。

不难发现 \(f(n)\) 其实是计算了 \(1 \sim n\) 每个数小于等于 \(n\) 的倍数个数之和，也就是说一个数对 \(f(n)\) 做出的贡献就是它的约数个数。那么有：

\[f(n)=\sum_{i=1}^{n}\lfloor\frac{n}{i}\rfloor=\sum_{i=1}^n{\rm{d}}(i) \]

由于 \({\rm{d}}\) 是积性函数，可以使用线性筛做到 \(O(n)\) 预处理。

关于线性筛 \({\rm{d}}\) ，可以看这个。

\({\rm{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\n",calcu(n,m));
	}
	return 0;
}