Solution

求
\[ \sum_{i=1}^{n}\sum_{j=1}^{m} lcm(i,j) \]
原式等价于
\[ \sum_{i=1}^{n}\sum_{j=1}^{m} \frac{ij}{gcd(i,j)} \]
老套路，令 \(gcd(i,j)=d\)，枚举一下 \(d\)
\[ \sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m} [gcd(i,j)=d] \times \frac{ij}{d} \]
毫无疑问地提出 \(d\)，\(d\) 对 \(\frac{ij}{d}\) 有贡献
\[ \sum_{d=1}^{n}d\times \sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}} [gcd(i,j)=1] \times ij \]
不用过脑子就把 \([gcd(i,j)=1]\) 换掉
\[ \sum_{d=1}^{n}d\times \sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}} \sum_{k|gcd(i,j)} \mu(k) \times ij \]
套路地枚举 \(k\)
\[ \sum_{d=1}^{n}d \times \sum_{k=1}^{\frac{n}{d}} \mu(k) \times \sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}} [k|gcd(i,j)] \times ij \]
继续
\[ \sum_{d=1}^{n}d \times \sum_{k=1}^{\frac{n}{d}} \mu(k) \times \sum_{ik=1}^{\frac{n}{d}}\sum_{jk=1}^{\frac{m}{d}}\times ij \times k^2 \]

\[ =\sum_{d=1}^{n}d \times \sum_{k=1}^{\frac{n}{d}} \mu(k) \times \sum_{i=1}^{\frac{n}{dk}}\sum_{j=1}^{\frac{m}{dk}}\times ij \times k^2 \]

归位
\[ =\sum_{d=1}^{n}d \times \sum_{k=1}^{\frac{n}{d}} \mu(k)k^2 \times \sum_{i=1}^{\frac{n}{dk}} \sum_{j=1}^{\frac{m}{dk}}\times ij \]
用等差数列替换
\[ =\sum_{d=1}^{n}d \times \sum_{k=1}^{\frac{n}{d}} \mu(k)k^2 \times \frac{\frac{n}{dk}(\frac{n}{dk} + 1)}{2} \times \frac{\frac{m}{dk}(\frac{m}{dk} + 1)}{2} \]
筛一下 \(\mu(k)k^2\)，然后数论分块套数论分块。

复杂度约为 \(O(n)\)

要到处膜。

Code

#include <bits/stdc++.h>
using namespace std;
#define re register
#define F first
#define S second
typedef long long ll;
typedef pair<int, int> P;
const int N = 1e7 + 5, M = 1e7 + 5, mod = 20101009;
const int INF = 0x3f3f3f3f;
inline int read() {
    int X = 0,w = 0; char ch = 0;
    while(!isdigit(ch)) {w |= ch == '-';ch = getchar();}
    while(isdigit(ch)) X = (X << 3) + (X << 1) + (ch ^ 48),ch = getchar();
    return w ? -X : X;
}
int p[N], mu[N];
ll sum[N];
bool vis[N];
int n, m;
void init(){
    int cnt = 0, k = min(n, m); mu[1] = 1;
    for (int i = 2; i <= k; i++){
        if (!vis[i]) p[++cnt] = i, mu[i] = -1;
        for (int j = 1; j <= cnt && i * p[j] <= k; j++){
            vis[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 <= N; i++) sum[i] = (sum[i - 1] + 1ll * i * i % mod * (mu[i] + mod) % mod) % mod;
}
ll get(int n, int m){
    return (1ll * n * (n + 1) / 2 % mod)  * (1ll * m * (m + 1) / 2 % mod) % mod;
}
ll solsum(int n, int m){
    ll ans = 0;
    for (int i = 1, j; i <= min(n, m); i = j + 1){
        j = min(n / (n / i), m / (m / i));
        ans = (ans + 1ll * (sum[j] - sum[i - 1] + mod) % mod * get(n / i, m / i) % mod) % mod;
    }
    return ans;
}
ll solve(int n, int m){
    ll ans = 0;
    for (int i = 1, j; i <= min(n, m); i = j + 1){
        j = min(n / (n / i), m / (m / i));
        ans = (ans + 1ll * (j - i + 1) * (i + j) / 2 % mod * solsum(n / i, m / i) % mod) % mod;
    }
    return ans;
}
int main(){
    n = read(), m = read();
    init();
    printf("%lld\n", solve(n, m));
    return 0;
}