求：
\(S(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}lcm(i,j)\)

显然：
\(S(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{ij}{gcd(i,j)}\)

枚举g：
\(S(n,m)=\sum\limits_{g=1}^{n}\frac{1}{g}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[gcd(i,j)==g]\)

除以g：
\(S(n,m)=\sum\limits_{g=1}^{n}g\sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{g}\rfloor}ij[gcd(i,j)==1]\)

记：
\(S_1(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[gcd(i,j)==1]\)

原式：
\(S(n,m)=\sum\limits_{g=1}^{n}gS_1(\lfloor\frac{n}{g}\rfloor,\lfloor\frac{m}{g}\rfloor)\)

化简\(S_1(n,m)\)，显然：
\(S_1(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij\sum\limits_{k|gcd(i,j)}\mu(k)\)

枚举k：
\(S_1(n,m)=\sum\limits_{k=1}^{min}\mu(k)\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[k|gcd(i,j)]\)

显然：
\(S_1(n,m)=\sum\limits_{k=1}^{min}\mu(k)\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[k|i][k|j]\)

这种时候可以除以k：
\(S_1(n,m)=\sum\limits_{k=1}^{min}\mu(k)k^2\sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{k}\rfloor}ij[1|i][1|j]\)

即：
\(S_1(n,m)=\sum\limits_{k=1}^{min}\mu(k)k^2\sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{k}\rfloor}ij\)

记：
\(S_2(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij\)

原式：
\(S_1(n,m)=\sum\limits_{k=1}^{min}\mu(k)k^2S_2(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor)\)

显然：
\(S_2(n,m)=\sum\limits_{i=1}^{n}i\sum\limits_{j=1}^{m}j\)

即：
\(S_2(n,m)=\frac{1}{4}n(n+1)m(m+1)\)

时间复杂度：
求\(S_2(n,m)\)是\(O(1)\)，分块求\(S_1(n,m)\)是\(O(n^{\frac{1}{2}})\)（大概），分块求\(S(n,m)\)是\(O(n)\)（大概）。
还需要线性筛出：\(\sum\limits_{k=1}^{min}\mu(k)k^2\)

线性筛已经足够了，还好写，不过为什么不试一波杜教筛呢？

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;const int mod=20101009;
const int MAXN=1e7;int pri[MAXN+1];
int &pritop=pri[0];
int A[MAXN+1];
int pk[MAXN+1];void sieve(int n=MAXN) {pk[1]=1;A[1]=1;for(int i=2; i<=n; i++) {if(!pri[i]) {pri[++pritop]=i;pk[i]=i;ll tmp=1ll*i*i;if(tmp>=mod)tmp%=mod;tmp=-tmp;if(tmp<mod)tmp+=mod;A[i]=tmp;}for(int j=1; j<=pritop; j++) {int &p=pri[j];int t=i*p;if(t>n)break;pri[t]=1;if(i%p) {pk[t]=p;ll tmp=1ll*A[i]*A[p];if(tmp>=mod)tmp%=mod;A[t]=tmp;} else {pk[t]=pk[i]*p;if(pk[t]==t) {A[t]=0;} else {ll tmp=1ll*A[t/pk[t]]*A[pk[t]];if(tmp>=mod)tmp%=mod;A[t]=tmp;}break;}}}for(int i=1; i<=n; i++) {A[i]+=A[i-1];if(A[i]>=mod)A[i]-=mod;}
}inline int qpow(ll x,int n) {ll res=1;while(n) {if(n&1) {res*=x;if(res>=mod)res%=mod;}x*=x;if(x>=mod)x%=mod;n>>=1;}return res;
}const int inv4=qpow(4,mod-2);inline int S2(int n,int m) {ll res=1ll*n*(n+1);if(res>=mod)res%=mod;res*=m;if(res>=mod)res%=mod;res*=(m+1);if(res>=mod)res%=mod;res*=inv4;if(res>=mod)res%=mod;return res;
}inline int S1(int n,int m) {ll res=0;int nm=min(n,m);for(int l=1,r; l<=nm; l=r+1) {int tn=n/l;int tm=m/l;r=min(n/tn,m/tm);ll tmp=A[r]-A[l-1];if(tmp<0)tmp+=mod;tmp*=S2(tn,tm);if(tmp>=mod)tmp%=mod;res+=tmp;}if(res>=mod)res%=mod;return res;
}inline int s1(int l,int r) {ll res=(1ll*(l+r)*(r-l+1))>>1;if(res>=mod)res%=mod;return res;
}inline int S(int n,int m) {ll res=0;int nm=min(n,m);for(int l=1,r; l<=nm; l=r+1) {int tn=n/l;int tm=m/l;r=min(n/tn,m/tm);ll tmp=s1(l,r);tmp*=S1(tn,tm);if(tmp>=mod)tmp%=mod;res+=tmp;}if(res>=mod)res%=mod;return res;
}int main() {
#ifdef Yinkufreopen("Yinku.in","r",stdin);
#endif // Yinkuint n,m;scanf("%d%d",&n,&m);sieve(max(n,m));printf("%d\n",S(n,m));return 0;
}

杜教筛还是非常快的。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;const int mod=20101009;
//杜教筛
const int MAXN=pow(1e7,2.0/3.0)+1;
int pri[MAXN+1];
int &pritop=pri[0];
int A[MAXN+1];
int pk[MAXN+1];void sieve(int n=MAXN) {pk[1]=1;A[1]=1;for(int i=2; i<=n; i++) {if(!pri[i]) {pri[++pritop]=i;pk[i]=i;ll tmp=1ll*i*i;if(tmp>=mod)tmp%=mod;tmp=-tmp;if(tmp<mod)tmp+=mod;A[i]=tmp;}for(int j=1; j<=pritop; j++) {int &p=pri[j];int t=i*p;if(t>n)break;pri[t]=1;if(i%p) {pk[t]=p;ll tmp=1ll*A[i]*A[p];if(tmp>=mod)tmp%=mod;A[t]=tmp;} else {pk[t]=pk[i]*p;if(pk[t]==t) {A[t]=0;} else {ll tmp=1ll*A[t/pk[t]]*A[pk[t]];if(tmp>=mod)tmp%=mod;A[t]=tmp;}break;}}}for(int i=1; i<=n; i++) {A[i]+=A[i-1];if(A[i]>=mod)A[i]-=mod;}
}inline int qpow(ll x,int n) {ll res=1;while(n) {if(n&1) {res*=x;if(res>=mod)res%=mod;}x*=x;if(x>=mod)x%=mod;n>>=1;}return res;
}const int inv4=qpow(4,mod-2);inline int S2(int n,int m) {ll res=1ll*n*(n+1);if(res>=mod)res%=mod;res*=m;if(res>=mod)res%=mod;res*=(m+1);if(res>=mod)res%=mod;res*=inv4;if(res>=mod)res%=mod;return res;
}const int inv6=qpow(6,mod-2);inline int s2(int n) {ll res=1ll*n*(n+1);if(res>=mod)res%=mod;res*=(2ll*n+1);if(res>=mod)res%=mod;res*=inv6;if(res>=mod)res%=mod;return res;
}unordered_map<int,int> GA;
inline int Get_A(int n){if(n<=MAXN)return A[n];if(GA.count(n))return GA[n];ll res=1;for(int l=2,r;l<=n;l=r+1){int t=n/l;r=n/t;ll tmp=s2(r)-s2(l-1);if(tmp<0)tmp+=mod;tmp*=Get_A(t);if(tmp>=mod)tmp%=mod;res-=tmp;}res%=mod;if(res<0)res+=mod;return GA[n]=res;
}inline int S1(int n,int m) {ll res=0;int nm=min(n,m);for(int l=1,r; l<=nm; l=r+1) {int tn=n/l;int tm=m/l;r=min(n/tn,m/tm);ll tmp=Get_A(r)-Get_A(l-1);if(tmp<0)tmp+=mod;tmp*=S2(tn,tm);if(tmp>=mod)tmp%=mod;res+=tmp;}if(res>=mod)res%=mod;return res;
}inline int s1(int l,int r) {ll res=(1ll*(l+r)*(r-l+1))>>1;if(res>=mod)res%=mod;return res;
}inline int S(int n,int m) {ll res=0;int nm=min(n,m);for(int l=1,r; l<=nm; l=r+1) {int tn=n/l;int tm=m/l;r=min(n/tn,m/tm);ll tmp=s1(l,r);tmp*=S1(tn,tm);if(tmp>=mod)tmp%=mod;res+=tmp;}if(res>=mod)res%=mod;return res;
}int main() {
#ifdef Yinkufreopen("Yinku.in","r",stdin);
#endif // Yinkuint n,m;scanf("%d%d",&n,&m);sieve();printf("%d\n",S(n,m));return 0;
}