Sol
首先有個結論
\(\sum_{i=1}^{m}\sum_{j=1}^{n}d(i*j)=\sum_{i=1}^{m}\sum_{j=1}^{n}\sum_{x|i}\sum_{y|i}[gcd(x,y)==1]\)
證明：可以看po姐的博客
接著這個式子推
\[ 原式=\sum_{x=1}^{n}\sum_{y=1}^{m}([gcd(x, y)==1] * \sum_{x|i}\sum_{y|i} 1)\\ =\sum_{x=1}^{n}\sum_{y=1}^{m}[gcd(x, y)==1\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor]\\ 設f(i)=\sum_{x=1}^{n}\sum_{y=1}^{m}[gcd(x, y)==i\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor]\\ 設g(i)=\sum_{x|d}f(d) \]
f(i)可以通過莫比烏斯反演求出
考慮求g(i)
\[ g(i)=\sum_{i|gcd(x,y)}\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor\\ =\sum_{i|x}\sum_{i|y}\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor\\ =\sum_{x=1}^{\lfloor\frac{n}{i}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{y}\rfloor}\lfloor\frac{n}{x*i}\rfloor\lfloor\frac{m}{y*i}\rfloor\\ 換個元=\sum_{i=1}^{x}\sum_{j=1}^{y}\lfloor\frac{x}{i}\rfloor\lfloor\frac{y}{j}\rfloor\\ \]
這個東西\(\sum_{i=1}^{x}\lfloor\frac{x}{i}\rfloor\)就是每個數的倍數的個數和的和，就是每個數約數的個數和的和預處理一下，前綴和一下就好，于是每個g(i)就可以O(1) 求。。。（約數的個數是積性函數，也可以線性篩）
數論分塊什么的就不說了
# include <bits/stdc++.h>
# define RG register
# define IL inline
# define Fill(a, b) memset(a, b, sizeof(a))
using namespace std;
typedef long long ll;
const int _(5e4 + 1);IL ll Read(){char c = '%'; ll x = 0, z = 1;for(; c > '9' || c < '0'; c = getchar()) if(c == '-') z = -1;for(; c >= '0' && c <= '9'; c = getchar()) x = x * 10 + c - '0';return x * z;
}int prime[_], mu[_], d[_], pred[_], num;
bool isprime[_];IL void Prepare(){isprime[1] = 1; mu[1] = d[1] = 1;for(RG int i = 2; i < _; ++i){if(!isprime[i]){ prime[++num] = i; mu[i] = -1; d[i] = 2; pred[i] = 1; }for(RG int j = 1; j <= num && i * prime[j] < _; ++j){isprime[i * prime[j]] = 1;if(i % prime[j]){ mu[i * prime[j]] = -mu[i]; d[i * prime[j]] = d[i] * 2; pred[i * prime[j]] = 1; }else{mu[i * prime[j]] = 0;pred[i * prime[j]] = pred[i] + 1;d[i * prime[j]] = d[i] / (pred[i] + 1) * (pred[i] + 2);break;}}d[i] += d[i - 1]; mu[i] += mu[i - 1];}
}int main(RG int argc, RG char *argv[]){Prepare();for(RG int T = Read(); T; --T){RG int n = Read(), m = Read(); RG ll ans = 0;if(n > m) swap(n, m);for(RG int i = 1, j; i <= n; i = j + 1){j = min(n / (n / i), m / (m / i));ans += 1LL * (mu[j] - mu[i - 1]) * d[n / i] * d[m / i];}printf("%lld\n", ans);}return 0;
}
