$$lcm(i,j)=\frac{ij}{gcd(i,j)}$$

$$\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_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d]\frac{ij}{d}$$

$$=\sum_{d=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=1]ijd$$

$$=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=1]ij$$

$$f(x)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=x]ij$$

$$g(x)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[x|gcd(i,j)]ij$$

$$g(x)=\sum_{x|d}f(d)$$

$$g(x)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[x|gcd(i,j)]ij$$

$$g(x)=\sum_{i=1}^{\lfloor\frac{n}{dx}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{dx}\rfloor}ij \cdot x^2$$

$$g(x)=x^2\sum_{i=1}^{\lfloor\frac{n}{dx}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{dx}\rfloor}ij$$

$$ans=\sum_{d=1}^{n}d\cdot f(1)$$

$$f(1)=\sum_{d=1}^{\lfloor\frac{n}{}\rfloor}\mu(d)g(d)$$

$$ans=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=1]ij$$

$$f(1)=\sum_{i=1}^{n}\mu(i)g(i)$$

$$f(1)=\sum_{i=1}^{n}\mu(i)i^2\sum_{a=1}^{\lfloor\frac{n}{di}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{di}\rfloor}ab$$

#### Code-$O(n)$

#include<cmath>
#include<cstdio>
#include<iostream>
#include<algorithm>

#define ll long long
#define MOD 20101009
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))

const int N=1e7+2;
const int inf=1e9+9;

int mui[N],sum[N];
int vis[N],prime[N],cnt;

inline void _pre_mui(int n){
mui[1]=1;
for(int i=2;i<=n;++i){
if(!vis[i])prime[++cnt]=i,mui[i]=-1;
for(int j=1;j<=cnt;++j){
if(i*prime[j]>n)break;
vis[i*prime[j]]=1;
if(!(i%prime[j])){mui[i*prime[j]]=0;break;}
mui[id*prime[j]]=-mui[i];
}
}return;
}

inline int solve(int n,int m){
ll ans=0;
for(int l=1,r=0;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
int res=1ll*(1ll*(n/l)*(n/l+1)/2%MOD)*(1ll*(m/l)*(m/l+1)/2%MOD)%MOD;
ans+=1ll*(sum[r]-sum[l-1])%MOD*res%MOD;
ans%=MOD;
}return (ans+MOD)%MOD;
}

int main(){
int n,m,ans=0;
scanf("%d%d",&n,&m);
if(n>m)n^=m^=n^=m;
_pre_mui(n);
for(int i=1;i<=n;++i)
sum[i]=(sum[i-1]+1ll*i*i%MOD*mui[i]%MOD+MOD)%MOD;
for(int l=1,r=0;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
int res=1ll*(l+r)*(r-l+1)/2%MOD;
ans=(ans+1ll*solve(n/l,m/l)*res%MOD)%MOD;
}
printf("%d\n",ans);
return 0;
}


“哇咔咔咔卡掉你们的 O(n) !”

$O(\sqrt{n})$ 是过的了的，故考虑向着这个方向前进。

$$ans=\sum_{d=1}^{n}d\sum_{i=1}^{n}\mu(i)i^2\sum_{a=1}^{\lfloor\frac{n}{di}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{di}\rfloor}ab$$

$$ans=\sum_{T=1}^{n}\sum_{a=1}^{\lfloor\frac{n}{T}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{T}\rfloor}ab\sum_{i|T}i^2\frac{T}{i}\mu(i)$$

？？？？？？？

$$\sum_{d=1}^{n}d\sum_{i=1}^{n}\mu(i)i^2$$

• $k$ 是质数，这下子后面的 $i$ 只能是 $1$ 和 $k$ ，$1$ 的时候的值就是 $k$ ，$k$ 的时候的值是 $k^2\cdot \frac{k}{k}\cdot \mu(k)$ ，很显然这个时候的 $\mu(k)$ 的值是 $-1$ ，于是这个时候的值是 $-k^2$ ，那么这个时候 $sum[k]$ 的值是 $k-k^2$ 。

• $\mu(k)$ 为 $0$ ，这下子的话就肯定有一个 $j$ ，使得 $k$ 可以整除 $j^2$ ，这个时候假设就只能整除 $j^2$ ，也就是说 $\mu(k/j)$ 的值非 $0$ 。那么我们看看，在 $sum$ 所计算的式子中，只有 $T$ 的因子对 $T$ 产生贡献。考虑 $k/j$ 到 $k$ 多了什么因子。这个时候多的因子有两类，一类是包含了 $j^2$ 的，一类是只包含了 $j$ 的。第二类的可以先不管，因为之前 $k/j$ 中有了一个 $j$ ，这类因子的贡献已经算过了。那么对于第一类因子，因为包含了 $j^2$ ，所以 $\mu$ 值为 $0$ ，对答案没有任何贡献。

那么这个时候对答案有贡献的还是 $k/j$ 的因子，乘上一个 $j$ 后没有更多的对答案造成贡献的因子。

但是我们发现上限 $T$ 变了，增大了 $j$ 倍，对于原来的每份贡献的值也增大了 $j$ 倍。由于没有其他的贡献，$k$ 的所有的贡献都来自 $k/j$ ，那么直接转移就好。

所以是 $sum[k]=sum[k/j]\cdot j$

• 对于剩下的情况，我们发现，这个可以直接转移了。当我们枚举 $k$ 的时候，考虑怎么用 $k$ 来转移 $k\cdot j$ ，这个时候 $j$ 是质数，并且 $k$ 中不包含 $j$ ，也就是说 $k$ 与 $j$ 互质。于是根据积性函数的性质，$sum[k]=sum[k/j]\cdot sum[j]$ 就好。

#### Code-$O(\sqrt{n})$

#include<cmath>
#include<cstdio>
#include<iostream>
#include<algorithm>

#define ll long long
#define MOD 20101009
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))

const int N=1e7+3;
const int inf=1e9+9;

int sum[N],vis[N],prime[N],cnt;

inline void _pre_mui_sum(){
vis[1]=sum[1]=1;
for(int i=2;i<=N;++i){
if(!vis[i])prime[++cnt]=i,sum[i]=(i-1ll*i*i%MOD+MOD)%MOD;
for(int j=1;j<=cnt;++j){
if(i*prime[j]>=N)break;
vis[i*prime[j]]=1;
if(!(i%prime[j])){sum[i*prime[j]]=1ll*sum[i]*prime[j]%MOD;break;}
else sum[i*prime[j]]=1ll*sum[i]*sum[prime[j]]%MOD;
}
}
for(int i=1;i<=N;++i)
sum[i]=(sum[i-1]+sum[i])%MOD;
}

int main(){
int n,m;
_pre_mui_sum();
scanf("%d%d",&n,&m);
if(n>m)std::swap(n,m);
long long ans=0;
for(int l=1,r=0;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
int res=(1ll*(1+n/l)*(n/l)/2%MOD)*(1ll*(1+m/l)*(m/l)/2%MOD)%MOD;
ans+=1ll*(sum[r]-sum[l-1]+MOD)%MOD*res%MOD;
ans%=MOD;
}
printf("%lld\n",(ans+MOD)%MOD);
return 0;
}


QAQ