易知原式等价于
$_math_inline$ \sum_{i=1}^n\sum_{j=1}^m\frac{i\cdot j}{\gcd(i,j)} $math_inline_$枚举最大公因数 $_math_inline$d$math_inline_$ ，显然两个数除以 $_math_inline$d$math_inline_$ 得到的数互质
$_math_inline$ \sum_{i=1}^n\sum_{j=1}^m\sum_{d\mid i,d\mid j,\gcd(\frac{i}{d},\frac{j}{d})=1}\frac{i\cdot j}{d} $math_inline_$非常经典的 $_math_inline$\gcd$math_inline_$ 式子的化法
$_math_inline$ \sum_{d=1}^n d\cdot\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1]\ i\cdot j $math_inline_$后半段式子中，出现了互质数对之积的和，为了让式子更简洁就把它拿出来单独计算。于是我们记
$_math_inline$ \text{sum}(n,m)=\sum_{i=1}^n\sum_{j=1}^m [\gcd(i,j)=1]\ i\cdot j $math_inline_$接下来对 $_math_inline$\text{sum}(n,m)$math_inline_$ 进行化简。首先枚举约数，并将 $_math_inline$[\gcd(i,j)=1]$math_inline_$ 表示为 $_math_inline$\varepsilon(\gcd(i,j))$math_inline_$
$_math_inline$ \sum_{d=1}^n\sum_{d\mid i}^n\sum_{d\mid j}^m\mu(d)\cdot i\cdot j $math_inline_$设 $_math_inline$i=i'\cdot d$math_inline_$ ， $_math_inline$j=j'\cdot d$math_inline_$ ，显然式子可以变为
$_math_inline$ \sum_{d=1}^n\mu(d)\cdot d^2\cdot\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}i\cdot j $math_inline_$观察上式，前半段可以预处理前缀和；后半段又是一个范围内数对之和，记
$_math_inline$ g(n,m)=\sum_{i=1}^n\sum_{j=1}^m i\cdot j=\frac{n\cdot(n+1)}{2}\times\frac{m\cdot(m+1)}{2} $math_inline_$可以 $_math_inline$\Theta(1)$math_inline_$ 求解
至此
$_math_inline$ \text{sum}(n,m)=\sum_{d=1}^n\mu(d)\cdot d^2\cdot g(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor) $math_inline_$我们可以 $_math_inline$\lfloor\frac{n}{\lfloor\frac{n}{d}\rfloor}\rfloor$math_inline_$ 数论分块求解 $_math_inline$\text{sum}(n,m)$math_inline_$ 函数。
在求出 $_math_inline$\text{sum}(n,m)$math_inline_$ 后，回到定义 $_math_inline$\text{sum}$math_inline_$ 的地方，可得原式为
$_math_inline$ \sum_{d=1}^n d\cdot\text{sum}(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor) $math_inline_$可见这又是一个可以数论分块求解的式子！
本题除了推式子比较复杂、代码细节较多之外，是一道很好的莫比乌斯反演练习题！（上述过程中，默认 $_math_inline$n\leqslant m$math_inline_$ ）
时间复杂度： $_math_inline$\Theta(n+m)$math_inline_$ （两次数论分块）
Code
1#include <bits/stdc++.h>
2
3using namespace std;
4
5const int N = 1e7;
6const int mod = 20101009;
7int n, m, mu[N + 5], p[N / 10 + 5], sum[N + 5];
8bool flg[N + 5];
9
10void init() {
11 mu[1] = 1;
12 int tot = 0, k = min(n, m);
13 for (int i = 2; i <= k; ++i) {
14 if (!flg[i]) p[++tot] = i, mu[i] = -1;
15 for (int j = 1; j <= tot && i * p[j] <= k; ++j) {
16 flg[i * p[j]] = 1;
17 if (i % p[j] == 0) {
18 mu[i * p[j]] = 0;
19 break;
20 }
21 mu[i * p[j]] = -mu[i];
22 }
23 }
24 for (int i = 1; i <= k; ++i)
25 sum[i] = (sum[i - 1] + 1LL * i * i % mod * (mu[i] + mod)) % mod;
26}
27
28int Sum(int x, int y) {
29 return (1LL * x * (x + 1) / 2 % mod) * (1LL * y * (y + 1) / 2 % mod) % mod;
30}
31
32int func(int x, int y) {
33 int res = 0;
34 for (int i = 1, j; i <= min(x, y); i = j + 1) {
35 j = min(x / (x / i), y / (y / i));
36 res = (res + 1LL * (sum[j] - sum[i - 1] + mod) * Sum(x / i, y / i) % mod) % mod;
37 }
38 return res;
39}
40
41int solve(int x, int y) {
42 int res = 0;
43 for (int i = 1, j; i <= min(x, y); i = j + 1) {
44 j = min(x / (x / i), y / (y / i));
45 res = (res + 1LL * (j - i + 1) * (i + j) / 2 % mod * func(x / i, y / i) % mod) % mod;
46 }
47 return res;
48}
49
50int main() {
51 cin >> n >> m;
52 init();
53 cout << solve(n, m) << endl;
54}
除另有声明外，本博客文章均采用 知识共享 (Creative Commons) 署名 4.0 国际许可协议 进行许可。转载请注明原作者与文章出处。