两次数论分块。
式子：
[sum_{i=1}^{n}sum_{j=1}^{m} lcm(i,j)\
= sum_{i=1}^{n}sum_{j=1}^{m} frac{ij}{gcd(i,j)}\
= sum_{i=1}^{n}sum_{j=1}^{m} sum_{d=gcd(i,j)}frac{ij}{d}\
= sum_{i=1}^{n}sum_{j=1}^{m} sum_{d=gcd(i,j)}dfrac{i}{d}frac{j}{d}\
= sum_{d=1}^{n} d sum_{i=1}^{n}sum_{j=1}^{m} frac id frac jd[gcd(i,j) = d]\
= sum_{d=1}^{n} dsum_{i=1}^{frac nd} sum_{j=1}^{frac md} ij[gcd(i,j) =1]\
*\
*\
f(n,m) = sum_{i=1}^{n} sum_{j=1}^{m}ij[gcd(i,j) = 1]\
= sum_{i=1}^{n} sum_{j=1}^{m} ij sum_{p | gcd(i,j) } mu(p)\
= sum_{i=1}^{n} sum_{j=1}^{m} ij sum_{p | i, p | j } mu(p)\
=sum_{p=1}^{n} mu(p) p^2 sum_{i=1}^{frac np}sum_{j=1}^{frac mp} ij\
=sum_{p=1}^{n} mu(p) p^2 g(frac np) g(frac mp)\
g(n) = frac{n(n + 1)}{2}\
*\
*\
sum_{d=1}^{n} dsum_{i=1}^{frac nd} sum_{j=1}^{frac md} ij[gcd(i,j) =1]\
= sum_{d=1}^{n} d f(frac nd , frac md)\
]
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 1e7;
const int mod = 20101009;
int p[N + 5], prime[N / 10], cnt;
int mu[N + 5], sum[N + 5];
ll g(ll n) { return (1ll * n * (n + 1) / 2) % mod; }
ll f(ll n, ll m){
ll ret = 0;
for(int l = 1,r; l <= n; l = r + 1){
r = min(min(n / (n/l), m / (m/l)), n);
ll v = 1ll * g(n/l) * g(m/l) % mod;
ret = (ret + 1ll * ((sum[r] - sum[l - 1] + mod) % mod) * v % mod) % mod;
}
return ret;
}
int n,m;
int main(){
p[0] = p[1] = 1; mu[1] = 1;
for(int i = 2; i <= N; ++ i){
if(!p[i]) { prime[++ cnt] = i; mu[i] = -1; }
for(int j = 1; j <= cnt && 1ll * prime[j] * i <= N; ++ j){
p[prime[j] * i] = 1;
if(i % prime[j] == 0) { mu[prime[j] * i] = 0; break; }
else mu[prime[j] * i] = - mu[i];
}
}
for(int i = 1; i <= N; ++ i) sum[i] = (sum[i - 1] + (1ll * mu[i] * (1ll * i * i % mod) % mod) % mod + mod) % mod;
scanf("%d%d",&n,&m);
if(n > m) swap(n,m);
ll ans = 0;
for(int l = 1, r; l <= n; l = r + 1){
r = min(min(n / (n/l), m / (m/l)), n);
ll v = f(n/l, m/l);
ans = (ans + 1ll * v * ((1ll * (r - l + 1) * (l + r) / 2) % mod) % mod ) % mod;
}
printf("%lld
",ans);
return 0;
}