大致题意： 求(sum_{i=1}^nsum_{j=1}^mlcm(i,j))。
推式子
不会莫比乌斯反演的可以先去看这篇博客：初学莫比乌斯反演。
反演题显然就是推式子啊~~~
考虑(lcm)这个东西不好做，所以我们先把它化成(gcd)：
[sum_{i=1}^nsum_{j=1}^mfrac{ij}{gcd(i,j)}
]
接下来就是按照套路枚举(gcd)了：
[sum_{d=1}^{min(n,m)}dsum_{i=1}^{lfloorfrac nd
floor}sum_{j=1}^{lfloorfrac md
floor}ij[gcd(i,j)=1]
]
由于([gcd(i,j)=1])这个式子不好直接搞，因此我们要进行转化。
想到(e(n)=[n=1])，所以这个式子就相当于(e(gcd(i,j)))。
而(mu*I=e)，所以：
[e(gcd(i,j))=sum_{p|gcd(i,j)}mu(p)cdot I(frac{gcd(i,j)}p)=sum_{p|gcd(i,j)}mu(p)
]
也就是说，原式即为：
[sum_{d=1}^{min(n,m)}dsum_{i=1}^{lfloorfrac nd
floor}sum_{j=1}^{lfloorfrac md
floor}ijsum_{p|gcd(i,j)}mu(p)
]
然后我们改为枚(p)，得到：
[sum_{d=1}^{min(n,m)}dsum_{p=1}^{lfloorfrac{min(n,m)}d
floor}sum_{i=1}^{lfloorfrac n{dp}
floor}sum_{j=1}^{lfloorfrac m{dp}
floor}ijp^2mu(p)
]
整理一下式子，调整一下顺序，得到：
[sum_{d=1}^{min(n,m)}dsum_{p=1}^{lfloorfrac{min(n,m)}d
floor}p^2mu(p)sum_{i=1}^{lfloorfrac n{dp}
floor}isum_{j=1}^{lfloorfrac m{dp}
floor}j
]
设(S(n)=sum_{i=1}^ni=frac{n(n+1)}2)，则可得原式等于：
[sum_{d=1}^{min(n,m)}dsum_{p=1}^{lfloorfrac{min(n,m)}d
floor}p^2mu(p)S(lfloorfrac n{dp}
floor)S(lfloorfrac m{dp}
floor)
]
那么我们预处理(p^2mu(p))的值，然后二维除法分块枚举(d,p)，就可以得出答案了。
代码
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 10000000
#define X 20101009
#define min(x,y) ((x)<(y)?(x):(y))
using namespace std;
int n,m,t,f[N+5];
class LinearSiever//线性筛
{
private:
int Pt,mu[N+5],P[N+5];
public:
int F[N+5];
I void Sieve(CI S)
{
RI i,j;for(mu[1]=1,i=2;i<=S;++i)//筛mu
{
!P[i]&&(mu[P[++Pt]=i]=-1);
for(j=1;j<=Pt&&1LL*i*P[j]<=S;++j)
if(P[i*P[j]]=1,i%P[j]) mu[i*P[j]]=-mu[i];else break;
}
for(i=1;i<=S;++i) F[i]=(1LL*i*i%X*(mu[i]+X)+F[i-1])%X;//预处理p^2mu(p)
}
}L;
I int S(CI x) {return (1LL*x*(x+1)>>1)%X;}//计算从1到x的和
int main()
{
RI l,r,ll,rr,res,ans=0;scanf("%d%d",&n,&m),L.Sieve(t=min(n,m));
for(l=1;l<=t;l=r+1)//第一维除法分块
{
res=0,r=min(n/(n/l),m/(m/l));
for(ll=1;ll<=t/l;ll=rr+1) rr=min((n/l)/(n/l/ll),(m/l)/(m/l/ll)),//第二维除法分块
res=(1LL*(L.F[rr]-L.F[ll-1]+X)%X*S(n/l/ll)%X*S(m/l/ll)+res)%X;
ans=(1LL*(S(r)-S(l-1)+X)%X*res+ans)%X;//统计答案
}return printf("%d",ans),0;//输出答案
}