这个题……主要问题在于(d(ij))应该怎么变形……容易想到改变成gcd的形式，不过不知道怎么改……
后来听大佬说有这么一个性质：
[d(ij) = sum_{x|i}sum_{y|j}[gcd(x,y) = 1]
]
这个不知道怎么严格证明……不过可以感性理解一下，就是首先肯定是要枚举i，j的因子的，然后如果这两个枚举的数不互质的话，那就可以把他们的gcd挪到一边使其变为互质的，也就是说如果二者不互质还计算，就会导致重复计数。（好吧这确实只是感性理解……）
那我们要求的就是
[sum_{i=1}^nsum_{j=1}^msum_{x|i}sum_{y|j}[gcd(x,y)=1]
]
这个式子的形式好像特别熟悉啊……把(sum)之间的限制互换，再用莫比乌斯函数的性质，再替换一下，这些熟悉的操作之后，这个式子就成了这样：
[sum_{d=1}^nmu(d)sum_{x=1}^{frac{n}{d}}sum_{y=1}^{frac{m}{d}}leftlfloorfrac{n}{dx}
ight
floorleftlfloorfrac{m}{dx}
ight
floor
]
然后发现外层的d可以用整除分块做，内层的可以先用整除分块预处理。然后就做完了。
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#include<set>
#include<vector>
#include<map>
#include<queue>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('
')
#define fr friend inline
#define y1 poj
#define mp make_pair
#define pr pair<int,int>
#define fi first
#define sc second
#define pb push_back
#define I puts("bug")
using namespace std;
typedef long long ll;
const int M = 50005;
const int INF = 1000000009;
const double eps = 1e-7;
int read()
{
int ans = 0,op = 1;char ch = getchar();
while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
return ans * op;
}
int n,m,T,p[M],mu[M];
ll pre[M],sum[M],tot;
bool np[M];
void euler()
{
np[1] = 1,mu[1] = 1;
rep(i,2,M-2)
{
if(!np[i]) p[++tot] = i,mu[i] = -1;
for(int j = 1;i * p[j] <= M-2 && j <= tot;j++)
{
np[i * p[j]] = 1;
if(!(i % p[j])) break;
mu[i * p[j]] = -mu[i];
}
}
rep(i,1,M-2) pre[i] = pre[i-1] + mu[i];
rep(k,1,M-2)
{
ll cur = 0;
for(int i = 1,j;i <= k;i = j + 1)
{
j = k / (k / i);
cur += (ll)(j - i + 1) * (ll)(k / i);
}
sum[k] = cur;
}
}
ll solve(ll a,ll b)
{
ll lim = min(a,b),cur = 0;
for(int i = 1,j;i <= lim;i = j + 1)
{
j = min(a / (a / i),b / (b / i));
cur += (pre[j] - pre[i-1]) * sum[a / i] * sum[b / i];
}
return cur;
}
int main()
{
euler();
T = read();
while(T--) n = read(),m = read(),printf("%lld
",solve(n,m));
return 0;
}