$$ \begin{aligned} \sum_{i=1}^n\sum_{j=1}^mlcm(i,j)&=\sum_{i=1}^n\sum_{j=1}^m\frac{ij}{(i,j)}\\\\ &=\sum_{d=1}^nd\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}ij[(i,j)=1]\\\\ &=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}ij\sum_{t|(i,j)}\mu(t)\\\\ &=\sum_{d=1}^nd\sum_{t=1}^{\lfloor\frac nd\rfloor}t^2\mu(t)\sum_{i=1}^{\lfloor\frac{\lfloor\frac nd\rfloor}{t}\rfloor}\sum_{j=1}^{\lfloor\frac{\lfloor\frac md\rfloor}{t}\rfloor}ij\\\\ &=\sum_{d=1}^nd\sum_{t=1}^{\lfloor\frac nd\rfloor}t^2\mu(t)\sum_{i=1}^{\lfloor\frac n{dt}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{dt}\rfloor}ij\\\\ &=\sum_{d=1}^nd\sum_{t=1}^{\lfloor\frac nd\rfloor}t^2\mu(t)\sum_{i=1}^{\lfloor\frac n{dt}\rfloor}i\sum_{j=1}^{\lfloor\frac{m}{dt}\rfloor}j\\\\ \end{aligned}\\\\ $$

后面一部分$\sum_{t=1}^{\lfloor\frac nd\rfloor}t^2\mu(t)\sum_{i=1}^{\lfloor\frac n{dt}\rfloor}i\sum_{j=1}^{\lfloor\frac{m}{dt}\rfloor}j$可以使用数论分块来求复杂度$O(\sqrt n)$

而前面的一部分也可以使用数论分块一共$O(\sqrt n)$ 块，因此总的时间复杂度是$O(n+m)$的.

// CREATED BY MASHIROYUKI

#include <bits/stdc++.h>
using namespace std;
#define IO ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
//#pragma GCC optimize("O3")
#define rep(i, x, n) for(int i = x; i <= n; i ++ )
#define downrep(i, x, n) for(int i = n; i >= x; i -- )

template<typename T> void debug(T x) {cout << "#debug: " << x << endl;}
template<typename T> void debug(T x, T y) {cout << "#debug: " << x << " " << y << endl;}
template<typename T>
void read(T &x){
    x = 0;
    char ch = getchar();
    while(!isdigit(ch)) ch = getchar();
    while(isdigit(ch)){
        x = (x << 1) + (x << 3) + (ch ^ 48);
        ch = getchar();
    }
}
typedef long long ll;
typedef double db;
typedef __int128 lll;
const int mod = 20101009;
const int N = 1e7 + 10;


ll pr[N],sum[N],mu[N],m,n;
bool vis[N];

void pre() {
  ll x=n,tot=0;
    mu[1]=1;
    rep(i,2,x) {
        if(!vis[i]) pr[++tot]=i,mu[i]=-1;
        for(int j=1;j<=tot&&(ll)i*pr[j]<=x;j++) {
            vis[i*pr[j]]=true;
            if(!(i%pr[j])) break;
            mu[i*pr[j]]=-mu[i];
        }
    }
    rep(i,1,x) sum[i]=(sum[i-1]+(ll)i*i%mod*(mu[i]+mod)%mod)%mod;
}

ll f(ll n,ll m) {
    ll l,r,res=0;
    for(l=1;l<=n;l=r+1) {
        ll v1=n/l,v2=m/l;
        r=(v1&&v2)?min(n,min(n/v1,m/v2)):n;
        res=((res+((sum[r]-sum[l-1])%mod+mod)%mod*lll(v1+1)*(v1)/2%mod*lll(v2+1)*v2/2)%mod+mod)%mod;
    }
    return (res%mod+mod)%mod;
}

void solve()
{
    read(n),read(m);
    if(n>m) swap(n,m);
    pre();
    ll l,r,res=0;
    for(l=1;l<=n;l=r+1) {
        ll v1=n/l,v2=m/l;
        r=(v1&&v2)?min(n,min(n/v1,m/v2)):n;
        res=(lll(l+r)*(r-l+1)/2%mod*f(v1,v2)+res)%mod;
    }
    printf("%lld\n",res);
}

signed main()
{
    //IO
    solve();
    return 0;
}