\(\frak{Solution}\)
\(\text{Update on 2022.2.6}\)：
发现自己写了好多废话，但是舍不得删，所以就变成 \(\rm details\) 了。
首先我们写出答案的式子（我们令 \(\mathtt{n <= m}\)）：
然后我们很自然地就想到了将 lcm 转化为 gcd。（毕竟公式是这么来的）
枚举 i 与 j 的 gcd。（P.S.\(\mathtt{[x=1]}\) 指的是当 x 取 1 时表达式的值为 1，其余为 0）
把 k 提出来。
我们讲一讲这个式子是怎么化过去的。我们知道，如果 \(\mathtt{gcd(i*k,j*k)==k}\)，那么其实可以转化成 \(\mathtt{gcd(i,j)==1}\)。可不可以这样理解：在 \(\mathtt{[1,n]}\) 与 \(\mathtt{[1,m]}\) 里分别选取 x，y，使其 gcd 为 k。在 \(\mathtt{[1,n/k]}\) 与 \(\mathtt{[1,m/k]}\) 中选择 i，j 使其互质。那么 \(\mathtt{i*k}\)，\(\mathtt{j*k}\) 和 \(\mathtt{x}\)，\(\mathtt{y}\) 其实是一一对应的关系。因为 \(\mathtt{n/k}\) 就是在 n 里面选取 k 的倍数，我们这样改写不会错过任何一个 gcd 为 k 的解。
至于后面的 \(\mathtt{i*j}\) 可以这样理解：通过前面的式子，我们知道后半部分其实就是 lcm。我们尝试列出此时的 lcm：\(\mathtt{\frac{(i*k)*(j*k)}{k}}\)。
即 \(\mathtt{i*j*k}\)，将 k 提前，就是 \(\mathtt{i*j}\)。
好了接下来我们利用莫比乌斯函数本身性质可以推出：
枚举 d，就是转化成 \(\mathtt{μ(d)}\) 可以匹配什么。。。至于为什么 d 的上界是 \(\mathtt{\frac nk}\)，我们将 \(\mathtt{gcd(i,j)}\) 都除了个 k 嘛。
我们枚举 i 或 j 倍的 d，把 \(\mathtt{i*k}\) 和 \(\mathtt{j*k}\) 看作原来的 i，j。因为肯定满足条件，所以直接乘起来。然后可以推出这样一个鬼式子（i，j 均为整）：
显然下一步可以化成
继续，
用等差数列求和公式：
最后我们筛出 \(\mathtt{μ(d)*d^2}\)，套两层数论分块就行了QwQ。
终于写完了。。。
emm 感觉证了半天自己对莫比乌斯反演一点感觉都没有我果然是蒟蒻。
首先对柿子进行转化：
考虑到莫反经常用到 "\([\gcd(i,j)=k]\)" 的形式，我们可以枚举 \(\gcd\)：
思考一下，如果令 \(g(k)=\sum_{i=1}^n\sum_{j=1}^m\ [\gcd(i,j)=k]\cdot \frac{i\times j}{k}\) 真的合理吗？不合理的原因在于后面带了个 \(\frac{1}{k}\)，由于 \(f(n)=\sum_{n\mid d}g(d)\)，此时显然不能快速求出 \(f\)（不确定是 \(\frac{1}{k}\) 还是 \(\frac{1}{3k}\) 之类……），也就达不到莫反优化的效果。
自然地，我们将 \(k\) 提出来：
这样可以观察到后面那一坨和 3.3. 用法叁 及其类似，直接套用可得：
枚举 \(t=d/k\)：
另外这里提醒一个需要 注意 的地方：一般莫反都是用 \(g(k)=\sum_{k\mid d}\mu(d/k)\cdot f(d)\)，为啥不写成 \(g(k)=\sum_{k\mid d}\mu(d)\cdot f(d/k)\)？
在这之前，我们先讨论 \(t\) 的取值范围从何而来：由于函数 \(g\) 的自变量取值 \(\le n\)，所以函数 \(f\) 的自变量取值 \(\le n\)，那么按照第一种莫反，函数 \(\mu\) 的自变量取值就 \(\le n/k\). 而对于第二种，函数 \(\mu\) 的自变量取值能达到惊人的 \(n^2\) 级别！
通常来看，我们知道 \(f,g\) 的自变量取值，所以本着减小复杂度的想法，一般使用第一种莫反。
再回到这道题，预处理前缀和，再用求和上的整除分块套上朴素整除分块即可。
\(\frak Code\)
#include <cstdio>
#include <cmath>
#include <iostream>
using namespace std;
#define rep(i,_l,_r) for(signed i=(_l),_end=(_r);i<=_end;++i)
#define print(x,y) write(x),putchar(y)
#define int long long
int read() {
int x=0,f=1; char s;
while((s=getchar())>'9'||s<'0') if(s=='-') f=-1;
while(s>='0'&&s<='9') x=(x<<1)+(x<<3)+(s^48),s=getchar();
return x*f;
}
void write(int x) {
if(x<0) return (void)(putchar('-'),write(-x));
if(x>9) write(x/10);
putchar(x%10^48);
}
const int maxn=1e7+5,mod=20101009;
int n,m,mu[maxn],p[(int)1e6],pc,ans,inv;
bool is[maxn];
int get(int x) {
return 1ll*(x+1)*x%mod*inv%mod;
}
int Duck(int N,int M) {
int r,ret=0,lim=min(N,M);
for(int l=1;l<=lim;l=r+1) {
r=min(N/(N/l),M/(M/l));
ret=(ret+1ll*(mu[r]-mu[l-1]+mod)%mod*get(N/l)%mod*get(M/l)%mod)%mod;
}
return ret;
}
void init() {
mu[1]=1;
rep(i,2,maxn-5) {
if(!is[i]) p[++pc]=i,mu[i]=-1;
rep(j,1,pc) {
if(p[j]*i>maxn-5) break;
is[p[j]*i]=1;
if(i%p[j]==0) break;
mu[p[j]*i]=-mu[i];
}
}
rep(i,1,maxn-5) mu[i]=(0ll+mu[i-1]+1ll*(mu[i]+mod)%mod*i%mod*i%mod)%mod;
}
int Goose() {
int lim=min(n,m),r,ret=0;
for(int l=1;l<=lim;l=r+1) {
r=min(n/(n/l),m/(m/l));
ret=(ret+1ll*(l+r)*(r-l+1)%mod*inv%mod*Duck(n/l,m/l)%mod)%mod;
}
return ret;
}
int qkpow(int x,int y) {
int r=1;
while(y) {
if(y&1) r=1ll*r*x%mod;
x=1ll*x*x%mod; y>>=1;
}
return r;
}
signed main() {
inv=qkpow(2,mod-2);
n=read(),m=read();
init();
printf("%lld\n",Goose());
return 0;
}
