题意

设\(d(x)\)为\(x\)的约数个数，求\(\sum_{i=1}^{n}\sum_{j=1}^{m}d(ij)\)。

题解

首先证个公式：

\[d(ij) = \sum_{x|i}\sum_{y|j} [gcd(x,y)=1]\]

可以这么考虑：利用唯一分解定理把\(i,j\)分解，即：

$i=\prod_{k = 1}^{m} p_k^{c_k},j=\prod_{k=1}^m p_k^{d_k} $

那等式左边显然为\(\prod(c_k+d_k+1)\)，

然后考虑等式右边在干什么事情：约数的最大公约数为\(1\)。

我们把\(x，y\)分解以后\(p_k\)的指数分别记为\(a_k,b_k\)

这就是说对于每个\(p_k\)，\(a_k,b_k\)最多有一个大于\(0\)。

那每一位分成三种情况，\(a_k=b_k=0\)；\(a_k\in[1,c_k]\)且\(b_k=0\)；\(a_k=0\)且\(b_k\in[1,d_k]\)

然后根据乘法原理乘起来，得到的结果和左边一样，就证出来了

然后带进去化简式子：（假设\(n\leq m\)）

\[=\sum_{i = 1}^n \sum_{j = 1}^m \sum_{x|i} \sum_{y|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 \]

\[=\sum_{x=1}^n \sum_{y=1}^m (\sum_{d|x,d|y} \mu(d) )\lfloor \frac{n}{x}\rfloor\lfloor \frac{m}{y} \rfloor\]

\[=\sum_{d=1}^n \mu(d) \sum_{x=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{y=1}^{\lfloor \frac{m}{d} \rfloor} \lfloor \frac{n}{xd} \rfloor \lfloor \frac{m}{xd} \rfloor\]

这里利用整除的小性质：\(\lfloor\frac{\lfloor \frac{a}{b} \rfloor}{c} \rfloor=\lfloor\frac{\lfloor\frac{a}{c}\rfloor}{b}\rfloor\)

\[=\sum_{d=1}^n \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\]

记\(f(n)=\sum_{i=1}^n\lfloor \frac{n}{i} \rfloor\)，则答案为\(\sum_{d=1}^n \mu(d)f(\lfloor \frac{n}{d} \rfloor) f(\lfloor \frac{m}{d} \rfloor)\)

这里如果\(f\)通过数论分块预处理，虽然可以过但是显然不够优美，可以换个角度考虑

这个式子的意义是，对于每个\(i\)统计它\(\leq n\)的所以倍数。那么每个数会被统计\(d(i)\)次（\(d\)为约数个数），因此：

\(f(n)=\sum_{i=1}^n d(i)\)

就可以线性筛预处理\(d\)的前缀和，还要预处理\(\mu\)的前缀和，然后每次数论分块。

#include <algorithm>
#include <cstdio>
using namespace std;typedef long long ll;const int N = 5e4 + 10;ll sum[N], mu[N], d[N];
int p[N], tot, e[N];
bool tag[N];void init(int n) {mu[1] = d[1] = 1;for(int i = 2; i <= n; i ++) {if(!tag[i]) {p[tot ++] = i; mu[i] = -1; d[i] = 2; e[i] = 1;}for(int j = 0; j < tot && i * p[j] <= n; j ++) {tag[i * p[j]] = 1;if(i % p[j] == 0) {mu[i * p[j]] = 0;e[i * p[j]] = e[i] + 1;d[i * p[j]] = d[i] / (e[i] + 1) * (e[i] + 2);break ;}mu[i * p[j]] = - mu[i];e[i * p[j]] = 1;d[i * p[j]] = d[i] * 2;}mu[i] += mu[i - 1];d[i] += d[i - 1];}
}ll calc(int n, int m) {ll ans = 0;for(int i = 1, j; i <= n; i = j + 1) {j = min(n / (n / i), m / (m / i));ans += (mu[j] - mu[i - 1]) * d[n / i] * d[m / i];}return ans;
}int main() {init(50000);int test; scanf("%d", &test);for(int n, m; test --; ) {scanf("%d%d", &n, &m);if(n > m) swap(n, m);printf("%lld\n", calc(n, m));}return 0;
}