$i \times j$的约数个数：

$$d(i, j) = \sum _ {x | i} \sum _ {y | i} [gcd(i, j) = 1]$$

\begin{aligned} & \sum _ {i = 1} ^ n \sum _ {j = 1} ^ m d(i, j)\\ =& \sum _ {i = 1} ^ n \sum _ {j = 1} ^ m \sum _ {x | i} \sum _ {y | i} [gcd(i, j) = 1]\\ =& \sum _ {i = 1} ^ n \sum _ {j = 1} ^ m \sum _ {x | i} \sum _ {y | i} \sum _ {t | x \text{&&} t | y} \mu(t)\\ =& \sum _ {t = 1} ^ n \sum _ {x = 1} ^ {\lfloor \frac n t \rfloor} \sum _ {y = 1} ^ {\lfloor \frac m t \rfloor} \sum _ {i = 1} ^ {\lfloor \frac n {tx} \rfloor} \sum _ y ^ {\lfloor \frac m {yt} \rfloor} \mu(t)\\ =& \sum _ {t = 1} ^ n \sum _ {x = 1} ^ {\lfloor \frac n t \rfloor} \sum _ {y = 1} ^ {\lfloor \frac m t \rfloor} \lfloor \frac n {tx} \rfloor \lfloor \frac m {yt} \rfloor \mu(t)\\ =& \sum _ {t = 1} ^ {n} (\sum _ {x = 1} ^ {\lfloor \frac n t \rfloor} {\lfloor \frac n {tx} \rfloor})(\sum _ {y = 1} ^ {\lfloor \frac m t \rfloor} {\lfloor \frac m {ty} \rfloor}) \end{aligned}

$$g(x) = \sum _ {i = 1} ^ x \lfloor \frac x i \rfloor$$

$$\sum _ {t = 1} ^ n g(\lfloor \frac n t \rfloor ) g(\lfloor \frac m t \rfloor )$$

（其实还有个前置知识：$\lfloor \frac a {bc} \rfloor = \lfloor \frac {\lfloor \frac a b \rfloor } c \rfloor$

#include <bits/stdc++.h>

#define NS (50005)
#define PS (6000)

using namespace std;

typedef long long LL;

template<typename _Tp> inline void IN(_Tp& dig)
{
char c; bool flag = 0; dig = 0;
while (c = getchar(), !isdigit(c)) if (c == '-') flag = 1;
while (isdigit(c)) dig = dig * 10 + c - '0', c = getchar();
if (flag) dig = -dig;
}

int testcase, n, m, p[PS], cnt, mu[NS], sum[NS];

bool vis[NS];

LL g[NS], ans;

void preWork()
{
mu[1] = 1;
for (int i = 2; i < NS; i += 1)
{
if (!vis[i]) p[++cnt] = i, mu[i] = -1;
for (int j = 1; j <= cnt; j += 1)
{
int k = i * p[j];
if (k >= NS) break;
vis[k] = 1, mu[k] = -mu[i];
if (i % p[j] == 0) { mu[k] = 0; break; }
}
}
for (int i = 1; i < NS; i += 1) sum[i] = sum[i - 1] + mu[i];
for (int i = 1; i < NS; i += 1)
for (int l = 1, r; l <= i; l = r + 1)
{
r = i / (i / l);
g[i] += 1ll * (r - l + 1) * (i / l);
}
}

int main(int argc, char const* argv[])
{
IN(testcase), preWork();
while (testcase--)
{
IN(n), IN(m), ans = 0;
if (n > m) swap(n, m);
for (int l = 1, r; l <= n; l = r + 1)
{
r = min(n / (n / l), m / (m / l));
ans += g[n / l] * g[m / l] * (sum[r] - sum[l - 1]);
}
printf("%lld\n", ans);
}
return 0;
}


#### Remmina

No puzzle that couldn't be solved.