洛谷

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

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

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

除以g：
(S(n,m)=sumlimits_{g=1}^{n}gsumlimits_{i=1}^{lfloorfrac{n}{g} floor}sumlimits_{j=1}^{lfloorfrac{m}{g} floor}ij[gcd(i,j)==1])

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

原式：
(S(n,m)=sumlimits_{g=1}^{n}gS_1(lfloorfrac{n}{g} floor,lfloorfrac{m}{g} floor))

化简(S_1(n,m))，显然：
(S_1(n,m)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}ijsumlimits_{k|gcd(i,j)}mu(k))

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

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

这种时候可以除以k：
(S_1(n,m)=sumlimits_{k=1}^{min}mu(k)k^2sumlimits_{i=1}^{lfloorfrac{n}{k} floor}sumlimits_{j=1}^{lfloorfrac{m}{k} floor}ij[1|i][1|j])

即：
(S_1(n,m)=sumlimits_{k=1}^{min}mu(k)k^2sumlimits_{i=1}^{lfloorfrac{n}{k} floor}sumlimits_{j=1}^{lfloorfrac{m}{k} floor}ij)

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

原式：
(S_1(n,m)=sumlimits_{k=1}^{min}mu(k)k^2S_2(lfloorfrac{n}{k} floor,lfloorfrac{m}{k} floor))

显然：
(S_2(n,m)=sumlimits_{i=1}^{n}isumlimits_{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))（大概）。
还需要线性筛出：(sumlimits_{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 Yinku
    freopen("Yinku.in","r",stdin);
#endif // Yinku
    int n,m;
    scanf("%d%d",&n,&m);
    sieve(max(n,m));
    printf("%d
",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 Yinku
    freopen("Yinku.in","r",stdin);
#endif // Yinku
    int n,m;
    scanf("%d%d",&n,&m);
    sieve();
    printf("%d
",S(n,m));
    return 0;
}