首先，把 $\operatorname{lcm}$ 转化成 $\gcd$ 可得原式等于 $\sum\limits_{i=1}^n\sum\limits_{j=1}^m\dfrac{ij}{\gcd(i,j)}$。

枚举 $d=\gcd(i,j)$ 可得上式等于 $\sum\limits_{d=1}^{\min(n,m)}d\sum\limits_{i=1}^{\left\lfloor\frac nd\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac md\right\rfloor}[\gcd(i,j)=1]\cdot ij$。

利用莫比乌斯反演，令 $f(n,m)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=1]\cdot ij=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{d|\gcd(i,j)}\mu(d)\cdot ij$。

调换枚举顺序，上式等于 $\sum\limits_{d=1}^{\min(n,m)}\mu(d)\sum\limits_{i=1}^{\left\lfloor\frac nd\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac md\right\rfloor}di\cdot dj=\sum\limits_{d=1}^{\min(n,m)}d^2\mu(d)\sum\limits_{i=1}^{\left\lfloor\frac nd\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac md\right\rfloor}ij$。后半部分可以根据 $\left\lfloor\dfrac nd\right\rfloor$ 和 $\left\lfloor\dfrac md\right\rfloor$ 的值直接计算，前半部分可以预处理出前缀和，所以整个式子的值可以用整除分块求解。

注意到我们要求 $\sum\limits_{d=1}^{\min(n,m)}d\cdot f\left(\left\lfloor\dfrac nd\right\rfloor,\left\lfloor\dfrac md\right\rfloor\right)$。显然，上式也可以整除分块求解。

总时间复杂度为 $O(\sqrt n\cdot\sqrt n)=O(n)$（假设 $n,m$ 同阶）。

/**
 * @date: 2024.01.22
 * @problem: P1829
 * @tags: 数学, 数论, 莫比乌斯反演, 整除分块
 */

#include<bits/stdc++.h>
using namespace std;
const int mod=20101009;

int n,m;

bool vis[10000001];
int mu[10000001],cntPrime,prime[1000001];
int sum[10000001];

void precalc(int n=10000000){
    mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!vis[i])prime[++cntPrime]=i,mu[i]=-1;
        for(int j=1;j<=cntPrime&&i*prime[j]<=n;j++){
            vis[i*prime[j]]=true;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;i++)
        sum[i]=(1ll*i*i%mod*mu[i]+sum[i-1]+mod)%mod;
}

inline int calcsum(int n,int m){
    return ((n+1ll)*n/2%mod)*((m+1ll)*m/2%mod)%mod;
}

int f(int n,int m){
    long long answer=0;
    for(int l=1,r;l<=n;l=r+1){
        r=min(n,min(n/(n/l),m/(m/l)));
        answer+=1ll*(sum[r]-sum[l-1])*calcsum(n/l,m/l)%mod;
    }
    return answer;
}

int main(){
    precalc();
    scanf("%d%d",&n,&m);
    if(n>m)swap(n,m);
    long long answer=0;
    for(int l=1,r;l<=n;l=r+1){
        r=min(n,min(n/(n/l),m/(m/l)));
        answer+=f(n/l,m/l)*((l+r)*(r-l+1ll)/2%mod)%mod;
    }
    printf("%lld\n",(answer%mod+mod)%mod);
    return 0;
}