提供两种（其实是一个式子）式子。
第一种
[sum_{g=1}^{min(n,m)} g sum_{d=1}^{min( lfloor frac{n}g
floor , lfloor frac{m}g
floor )} mu(d) * d^2 Bigg( sum_{i=1}^{lfloor frac{n}{dg}
floor} i sum_{j=1}^{lfloor frac{m}{dg}
floor} j Bigg)
]
第二种
[sum_{T=1}^{min(n,m)} Bigg( sum_{i=1}^{lfloor frac{n}{T}
floor} i sum_{j=1}^{lfloor frac{m}{T}
floor} j Bigg) sum_{d|T} d* mu(d) *T
]
其实都可以算， 时间所迫（还有很多题没做啊啊啊啊）， 只写第二种， 第一种口胡下 : 对 (g) 整除分块， 再对 (d) 整除分块， 复杂度是 (O(n^{frac{3}4}))， 再加上 (O(n)) 预处理 (mu)， 总复杂度是 (O(n)) 的。
第二种 (O(n)) 预处理出 (sum_{d|T} d* mu(d) *T) 的前缀和， 再分块 (O(sqrt n)) 算， 总复杂度也是 (O(n)) 的。
第二种Luogu数据AC代码：
#include<bits/stdc++.h>
using namespace std;
const int maxn = 1e7 + 15;
const int mod = 20101009;
const int inv2 = 10050505;
#define li long long
int pr_cnt, prime[maxn], v[maxn];
li f[maxn], s[maxn];
void euler(int n) {
f[1] = s[1] = 1ll;
for(int i=2; i<=n; ++i) {
if(!v[i]) {
v[prime[++pr_cnt] = i] = i;
f[i] = 1-i+mod; //不用取模qwq
}
for(int j=1; j<=pr_cnt; ++j) {
if(prime[j] > n/i || prime[j] > v[i]) break;
v[prime[j] * i] = prime[j];
f[prime[j] * i] = f[i] * (i%prime[j] ? f[prime[j]] : 1ll) % mod;
}
}
for(int i=1; i<=n; ++i) {
f[i] = (1ll*f[i]*i) % mod;
}
for(int i=2; i<=n; ++i) {
s[i] = s[i-1] + 1ll*i;
s[i] %= mod;
}
for(int i=1; i<=n; ++i) {
f[i] = (f[i]+f[i-1])%mod;
}
}
int main()
{
int n,m; cin>>n>>m; euler(max(n,m));
int len = min(n,m);
li ans = 0ll;
for(int i=1,j;i<=len;i=j+1) {
j = min(n/(n/i), m/(m/i));
j = min(j,len);
ans += s[n/i] * s[m/i] % mod * (f[j]-f[i-1]+mod) % mod;
ans = (ans%mod+mod)%mod;
}
cout << ans;
return 0;
}