要求的式子
$$\sum_{i = 1}^{N} \sum_{j = 1}^{M}d(i, j)$$
第一部将原式子转换
$$ \sum_{i = 1}^{N} \sum_{j = 1}^{M}\sum_{b|i}\sum_{d|j}[(b, d) = 1] $$
(上述的转换如果不懂没有办法，因为数论经常要用到，如果不会证明的可以去看y总的视频，或者直接百度)

要求什么就设什么为$f[]$
直接令
$$ f[x] =\sum_{i = 1}^{N} \sum_{j = 1}^{M}\sum_{b|i}\sum_{d|j}[(b, d) = x] $$
则 $f[1]$ 就是答案

根据反演公式 设出 F 考虑 F 的实际意义以及简便求法
$$ F[x] = \sum_{x|d}f[d] $$

$$ F[x] =\sum_{i = 1}^{N} \sum_{j = 1}^{M}\sum_{b|i}\sum_{d|j}[x|(b, d)] $$

发现 其实就是枚举 的b,d 都有 x 这个因子即可 考虑让 b, d 都除以 x

则有
$$ F[x] =\sum_{i = 1}^{N} \sum_{j = 1}^{M}\sum_{p|(i/x)}\sum_{q|(j/x)}1(i,j能够整除 x) $$
可以设
$$ D[x] =\sum_{q|x}1 $$
则
$$ F[x] =\sum_{i = 1}^{N} \sum_{j = 1}^{M}D[i/x]D[j/x] (这里有个条件, i,j能够整除 x) $$
转换枚举位置

$$ F[x] =\sum_{i = 1}^{N}D[i/x] \sum_{j = 1}^{M}D[j/x] $$
因为枚举过程中i, j 一定能够整除 x, 所以可以枚举倍数
$$ F[x] = \sum_{p = 1}^{N / x}D[p] \sum_{q = 1}^{M / x}D[q] $$
发现 可以通过预处理 D 函数 以及 D 的前缀和函数 那么可以O(1) 的时间求出 F 函数

然后F函数求出来之后再反演
$$ f[x] = \sum_{x|d}\mu[d/x]F[d] $$

$$ f[1] = \sum_{d}\mu[d]F[d] $$

$$ =\sum_{d}\mu[d]\sum_{p = 1}^{N / d}D[p] \sum_{q = 1}^{M / d}D[q] $$

发现又可以分块 那么就可以根号的时间出来
$$ 预处理时间O(n*\sqrt{n}) $$

$$ 查询时间O(t*\sqrt{t}) $$

$$ 总时间时间O((n + t)\sqrt{max(n + t)}) = O(n\sqrt{n}) $$

#include <bits/stdc++.h>
#define rep(i, a, b) for (int i = a; i <= b; ++i)
#define per(i, a, b) for (int i = a; i >= b; --i)
using namespace std;
typedef long long ll;
const int maxn = 5e4 + 5;
bool st[maxn];
int prime[maxn], cnt, mu[maxn], sm[maxn], D[maxn], Dm[maxn];
void init() {
    mu[1] = 1;
    rep(i, 2, maxn - 5) {
        if (!st[i]) prime[++cnt] = i, mu[i] = -1;
        rep(j, 1, cnt) {
            if ((ll)prime[j] * i > maxn - 5) break;
            st[prime[j] * i] = true;
            if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i];
        }
    }
    rep(i, 1, maxn - 5) {
        sm[i] = sm[i - 1] + mu[i];
        int tot = 0;
        rep(j, 1, i) {
            if ((ll)j * j > i) break;
            if (i % j == 0) {
                if (j == i / j) tot += 1;
                else tot += 2;
            }
        } 
        D[i] = tot;
    }
    rep(i, 1, maxn - 5) {
        Dm[i] = Dm[i - 1] + D[i];
    }
}
ll get(int n, int m) {
    ll ans = 0;
    for (int l = 1, r; l <= min(n, m); l = r + 1) {
        r = min(min(n, m), min(n / (n / l), m / (m / l)));
        ans += (ll)(sm[r] - sm[l - 1]) * Dm[n / l] * Dm[m / l];
    }
    return ans;
}
int main() {
    init();
    int t; cin >> t;
    while (t--) {
        int l, r;
        scanf("%d%d", &l, &r);
        printf("%lld\n", get(l, r));
    }
    return 0;
}