易知原式等价于
$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
#include <bits/stdc++.h>
using namespace std;
const int N = 1e7;
const int mod = 20101009;
int n, m, mu[N + 5], p[N / 10 + 5], sum[N + 5];
bool flg[N + 5];
void init() {
mu[1] = 1;
int tot = 0, k = min(n, m);
for (int i = 2; i <= k; ++i) {
if (!flg[i]) p[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot && i * p[j] <= k; ++j) {
flg[i * p[j]] = 1;
if (i % p[j] == 0) {
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
for (int i = 1; i <= k; ++i)
sum[i] = (sum[i - 1] + 1LL * i * i % mod * (mu[i] + mod)) % mod;
}
int Sum(int x, int y) {
return (1LL * x * (x + 1) / 2 % mod) * (1LL * y * (y + 1) / 2 % mod) % mod;
}
int func(int x, int y) {
int res = 0;
for (int i = 1, j; i <= min(x, y); i = j + 1) {
j = min(x / (x / i), y / (y / i));
res = (res + 1LL * (sum[j] - sum[i - 1] + mod) * Sum(x / i, y / i) % mod) % mod;
}
return res;
}
int solve(int x, int y) {
int res = 0;
for (int i = 1, j; i <= min(x, y); i = j + 1) {
j = min(x / (x / i), y / (y / i));
res = (res + 1LL * (j - i + 1) * (i + j) / 2 % mod * func(x / i, y / i) % mod) % mod;
}
return res;
}
int main() {
cin >> n >> m;
init();
cout << solve(n, m) << endl;
}