算法
(莫比乌斯反演) $O(t\sqrt{n})$
首先有一个结论:
$$ d(ij) = \sum_{x \mid i} \sum_{y \mid j} [\gcd(x, y) = 1] $$
证明这个结论有很多方式。注意到各个素因子的贡献是独立的。
而当 $i,j$ 都是某个素数的幂的时候这是显然的。
而对这个等式进行二维前缀和,我们有:
$$ \begin{aligned} \sum_{i = 1}^n \sum_{j = 1}^m d(ij) &= \sum_{i = 1}^n \sum_{j = 1}^m \sum_{x \mid i} \sum_{y \mid j} [\gcd(x,y) = 1] \\\ &= \sum_{x = 1}^n \sum_{y = 1}^m [\gcd(x,y) = 1] \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor \end{aligned} $$
利用 $\varepsilon = \mu * 1$,我们有
$$ \begin{aligned} \sum_{x = 1}^n \sum_{y = 1}^m [\gcd(x,y) = 1] \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor &= \sum_{x = 1}^n \sum_{y = 1}^m \left(\sum_{d \mid x, d \mid y} \mu(d)\right) \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor \\\ &= \sum_{d = 1}^{\min(n,m)} \mu(d) \sum_{x’ = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{y’ = 1}^{\lfloor \frac{m}{d} \rfloor} \lfloor \frac{n}{x’d} \rfloor \lfloor \frac{m}{y’d} \rfloor \\\ &= \sum_{d = 1}^{\min(n,m)} \mu(d) \sum_{x’ = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{y’ = 1}^{\lfloor \frac{m}{d} \rfloor} \lfloor \frac{\lfloor \frac{n}{d} \rfloor}{x’} \rfloor \lfloor \frac{\lfloor \frac{m}{d} \rfloor}{y’} \rfloor \end{aligned} $$
令 $\displaystyle f(n) = \sum_{i = 1}^n \lfloor \frac{n}{i} \rfloor$,则答案最终等于
$$ \sum_{d = 1}^{\min(n,m)} \mu(d) f(\lfloor \frac{n}{d} \rfloor) f(\lfloor \frac{m}{d} \rfloor) $$
而又有
$$ f(n) = \sum_{i = 1}^n \lfloor \frac{n}{i} \rfloor = \sum_{k = 1}^n d(k) $$
故可以利用 $\rm Euler$ 筛线性预处理 $d(n)$,从而前缀和得到 $f(n)$,之后就可以每次询问 $O(\sqrt{n})$ 回答了。
C++ 代码
#include <bits/stdc++.h>
using std::cin;
using std::cout;
using std::min;
using ll = long long;
const int N = 50010;
int primes[N], cnt, mu[N], sum[N], c[N], d[N], f[N];
bool st[N];
int g(int k, int x) {
return k / (k / x);
}
void init() {
mu[1] = 1, d[1] = 1;
for (int i = 2; i < N; ++i) {
if (!st[i]) primes[cnt++] = i, mu[i] = -1, c[i] = 1, d[i] = 2;
for (int j = 0; primes[j] * i < N; ++j) {
st[primes[j] * i] = true;
if (i % primes[j] == 0) {
d[i * primes[j]] = d[i] / (c[i] + 1) * (c[i] + 2);
c[i * primes[j]] = c[i] + 1;
break;
}
mu[i * primes[j]] = -mu[i];
d[i * primes[j]] = d[i] * d[primes[j]];
c[i * primes[j]] = 1;
}
}
for (int i = 1; i < N; ++i) sum[i] = sum[i - 1] + mu[i];
for (int i = 1; i < N; ++i) f[i] = f[i - 1] + d[i];
}
int main() {
init();
int t;
cin >> t;
while (t--) {
int n, m;
cin >> n >> m;
ll res = 0;
int k = min(n, m);
for (int l = 1, r; l <= k; l = r + 1) {
r = min(k, min(g(n, l), g(m, l)));
res += (ll)(sum[r] - sum[l - 1]) * f[n / l] * f[m / l];
}
cout << res << '\n';
}
return 0;
}