设无限集 \(\mathbb P\) 表示全体质数，则问题即求
\[\sum_{p \in \mathbb P}\sum_{x=1}^n \sum_{y=1}^m [\gcd(x,y)=p]
\]
发现有一个叫 \([\gcd(x,y)=p]\) 的东西，那不就是莫比乌斯反演吗？？？不会的话建议重修 P3327。
令 \(x'=\lfloor x/p \rfloor,\lfloor y'=y/p \rfloor\)，则原式等于
\[\sum_{p \in \mathbb P} \sum_{x'=1}^{\lfloor n/p \rfloor} \sum_{y'=1}^{\lfloor m/p \rfloor} [\gcd(x',y')=1]
\]
又是 \([k=1]\) 的形式，直接替换为莫比乌斯函数，即
\[\sum_{p \in \mathbb P} \sum_{x'=1}^{\lfloor x/p \rfloor} \sum_{y'=1}^{\lfloor m/p \rfloor} \sum_{d \mid \gcd(x',y')} \mu(d)
\]
交换求和顺序，先枚举 \(d\)，则有
\[\sum_{p \in \mathbb P} \sum_{d=1}^{\min(\lfloor x/p \rfloor,\lfloor y/p \rfloor)} \mu(d) \cdot \lfloor \frac n{pd} \rfloor \cdot \lfloor \frac m{pd} \rfloor
\]
我们发现有 \(pd\) 这个部分出现在了两个向下取整里，我们考虑交换求和顺序，先枚举 \(pd\)，即
\[\sum_{t=1}^{\min(n,m)} \lfloor \frac nt \rfloor \lfloor \frac mt \rfloor \sum_{p \mid t,p \in \mathbb P} \mu(t/p)
\]
记
\[f(t)=\sum_{p \mid t,p \in \mathbb P} \mu(t/p)
\]
答案可化为
\[\sum_{t=1}^{\min(n,m)} \lfloor \frac nt \rfloor \lfloor \frac mt \rfloor f(t)
\]
预处理 \(f(t)\) 的前缀和，剩下两个向下取整进行整除分块即可。
code
#include <bits/stdc++.h>
#define fast std::ios::sync_with_stdio(false); std::cin.tie(NULL); std::cout.tie(NULL);
constexpr std::int32_t mod = 998244353;
constexpr std::int32_t maxn = 1e7 + 5;
std::int32_t n, m, tot;
std::int64_t ans;
std::int32_t prime[maxn], mu[maxn];
std::int64_t f[maxn], pre_fsum[maxn];
bool vis[maxn];
int main() {std::freopen("contest.in", "r", stdin);std::freopen("contest.out", "w", stdout);fast;mu[1] = 1;for (std::int32_t i = 2; i <= maxn - 5; i++){if (!vis[i]) prime[++tot] = i, mu[i] = -1, f[i] = 1;for (std::int32_t j = 1; j <= tot && i * prime[j] <= maxn - 5; j++){vis[i * prime[j]] = true;if (i % prime[j] == 0){mu[i * prime[j]] = 0, f[i * prime[j]] = mu[i];break;} else mu[i * prime[j]] = -mu[i], f[i * prime[j]] = mu[i] - f[i];}}for (std::int32_t i = 1; i <= maxn - 5; i++) pre_fsum[i] = pre_fsum[i - 1] + f[i];std::int32_t t; std::cin >> t;while (t--){std::cin >> n >> m;if (n > m) std::swap(n, m); ans = 0;for (std::int32_t pl = 1, pr; pl <= n; pl = pr + 1){pr = std::min(n / (n / pl), m / (m / pl));ans += (std::int64_t)(n / pl) * (m / pl) * (pre_fsum[pr] - pre_fsum[pl - 1]);}std::cout << ans << std::endl;}return 0;
}