算法
(莫比乌斯反演) $O(t\sqrt{n})$
我们先考虑:使得 $\gcd(x, y) = 1$,即互质的 $x,y$ 组数有多少?
$$ \begin{aligned} \sum_{x = 1}^n \sum_{y = 1}^m [\gcd(x, y) = 1] &= \sum_{x = 1}^n \sum_{y = 1}^m \sum_{d \mid \gcd(x, y)} \mu(d) \\\ &= \sum_{x = 1}^n \sum_{y = 1}^m \sum_{d \mid x, d \mid y} \mu(d) \\\ &= \sum_{d = 1}^{\min(n, m)} \mu(d)\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor \end{aligned} $$
预处理 $\mu$ 及其前缀和,按 $\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor$ 的取值分段,就可以 $O(\sqrt{n})$ 解决这个问题。
记素数的集合为 $P$,那么原题就是又加了一重求和。
$$ \begin{aligned} \sum_{p \in P} \sum_{x = 1}^n \sum_{y = 1}^m [\gcd(x, y) = p] &= \sum_{p \in P} \sum_{x = 1}^{\lfloor \frac{n}{p} \rfloor} \sum_{y = 1}^{\lfloor \frac{m}{p} \rfloor} [\gcd(x, y) = 1] \\\ &= \sum_{p \in P} \sum_{x = 1}^{\lfloor \frac{n}{p} \rfloor} \sum_{y = 1}^{\lfloor \frac{m}{p} \rfloor} \sum_{d \mid x, d \mid y} \mu(d) \\\ &= \sum_{p \in P} \sum_{d = 1}^{\lfloor \frac{\min(n, m)}{p} \rfloor} \mu(d) \lfloor \frac{n}{pd}\rfloor \lfloor \frac{m}{pd} \rfloor \\\ &= \sum_{t = 1}^{\min(n, m)} \lfloor \frac{n}{t} \rfloor \lfloor \frac{m}{t} \rfloor \sum_{p \in P, p \mid t} \mu(\frac{t}{p}) \end{aligned} $$
令 $f(n) = [n \in P]$,并且
$$ g(n) = \sum_{p \in P, p \mid n} \mu(\frac{n}{p}) $$
那么有 $g = f * \mu$。
可以直接计算 $g$ 的前 $n$ 项及前缀和,之后就可以每次询问 $O(\sqrt{n})$ 回答了。
注意到计算 $g$ 的前 $n$ 项只需要枚举素数的倍数,这样做的时间复杂度实际上是 $O(n\log\log n)$。
C++ 代码
#include <bits/stdc++.h>
using std::cin;
using std::cout;
using ll = long long;
const int N = 10000010;
int primes[N], cnt;
bool st[N];
int g[N], u[N], sum[N];
void init() {
std::memset(st, 0, sizeof st);
u[1] = 1;
cnt = 0;
for (int i = 2; i < N; ++i) {
if (!st[i]) {
primes[cnt++] = i;
u[i] = -1;
g[i] = 1;
}
for (int j = 0; j < cnt and i * primes[j] < N; ++j) {
st[i * primes[j]] = 1;
if (i % primes[j]) {
u[i * primes[j]] = -u[i];
g[i * primes[j]] = u[i] - g[i];
}
else {
u[i * primes[j]] = 0;
g[i * primes[j]] = u[i];
break;
}
}
}
sum[0] = 0;
for (int i = 1; i < N; ++i) {
sum[i] = sum[i - 1] + g[i];
}
}
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
init();
int t;
cin >> t;
while (t--) {
ll n, m;
cin >> n >> m;
if (n > m) std::swap(n, m);
ll ans = 0;
for (int i = 1, last; i <= n; i = last + 1) {
last = std::min(n / (n / i), m / (m / i));
ans += (n / i) * (m / i) * (sum[last] - sum[i - 1]);
}
cout << ans << '\n';
}
return 0;
}