這道題我們要求的是

\[\sum_{i=1}^N\sum_{j=1}^Mlcm(i,j)\]

總所周知\(lcm\)的性質不如\(gcd\)優雅，但是唯一分解定理告訴我們\(gcd(i,j)\times lcm(i,j)=i\times j\)

所以很容易的可以轉化成這個柿子

\[\sum_{i=1}^N\sum_{j=1}^M\frac{i\times j}{(i,j)}\]

現在開始套路了

先設兩個函式

\[f(n)=\sum_{i=1}^N\sum_{j=1}^M[(i,j)==n]\times i\times j\]

\[F(n)=\sum_{i=1}^N\sum_{j=1}^M[n|(i,j)]\times i\times j\]

\[=\frac{n\times(\left \lfloor \frac{N}{n}\right \rfloor+1)\times \left \lfloor \frac{N}{n} \right \rfloor}{2}\times\frac{n\times(\left \lfloor \frac{M}{n}\right \rfloor+1)\times \left \lfloor \frac{M}{n} \right \rfloor}{2}\]

顯然則有

\[F(n)=\sum_{n|d}f(d)\]

反演得

\[f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d)\]

於是答案就是

\[Ans=\sum_{i=1}^N\frac{f(i)}{i}\]

\[=\sum_{i=1}^N\times\frac{1}{i}\sum_{i|d}\mu(\frac{d}{i})\frac{d\times(\left \lfloor \frac{N}{d}\right \rfloor+1)\times \left \lfloor \frac{N}{d} \right \rfloor}{2}\times\frac{d\times(\left \lfloor \frac{M}{d}\right \rfloor+1)\times \left \lfloor \frac{M}{d} \right \rfloor}{2}\]

後面的一大坨東西真是太煩人了，搞到前面來

\[Ans=\sum_{d=1}^N\frac{(\left \lfloor \frac{N}{d}\right \rfloor+1)\times \left \lfloor \frac{N}{d} \right \rfloor\times(\left \lfloor \frac{M}{d}\right \rfloor+1)\times \left \lfloor \frac{M}{d} \right \rfloor}{4}\sum_{i|d}\frac{\mu(\frac{d}{i})\times d^2}{i}\]

於是我們可以用\(\Theta(n \ln n)\)來求出\(\sum_{i|d}\frac{\mu(\frac{d}{i})\times d^2}{i}\)之後字首和

於是有了一個\(70\)分的程式碼

#include<iostream>
#include<cstring>
#include<cstdio>
#define re register
#define maxn 10000005
#define LL long long
#define min(a,b) ((a)<(b)?(a):(b))
#define max(a,b) ((a)>(b)?(a):(b))
const LL mod=20101009;
int n,m;
int f[maxn],p[maxn>>2];
LL pre[maxn];
int main()
{
    scanf("%d%d",&n,&m);
    if(n>m) std::swap(n,m);
    f[1]=pre[1]=1;
    for(re int i=2;i<=n;i++)
    {
        if(!f[i]) p[++p[0]]=i,pre[i]=(1-i+mod);
        for(re int j=1;j<=p[0]&&p[j]*i<=n;j++)
        {
            f[p[j]*i]=1;
            if(i%p[j]==0) 
            {
                pre[i*p[j]]=pre[i];
                break;
            }
            pre[p[j]*i]=pre[p[j]]*pre[i]%mod;
        }
    }
    for(re int i=1;i<=n;i++) pre[i]=(i*pre[i]%mod+pre[i-1])%mod;
    LL ans=0;
    for(re LL l=1,r;l<=n;l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        LL N=n/l,M=m/l;
        ans=(ans+((N+1)*N/2%mod)*((M+1)*M/2%mod)%mod*(pre[r]-pre[l-1]+mod)%mod)%mod;
    }
    printf("%lld\n",ans);
    return 0;
}

顯然不行

我們設\(T=\frac{d}{i}\)

那麼

\[\sum_{i|d}\frac{\mu(\frac{d}{i})\times d^2}{i}=d\times\sum_{T|d}\mu(T)\times T\]

定義\(h(x)=\sum_{T|x}\mu(T)\times T\)

會發現\(h\)是一個積性函式，可以考慮如何線篩

首先\(x\)是質數\(h(x)=1-x\)

互質的話可以直接乘起來

如果不互質的話需要好好考慮一下了

仔細思考一下這個\(h\)的含義，會發現有一些約數\(T\)是沒有用的，就是那些\(\mu(T)=0\)的約數

而線篩的時候一旦\(i\%p[j]==0\)，說明\(p[j]\)在\(i\)中出現過，於是\(p[j]\)並不能組成一些新的有用約數，函式值和\(h(i)\)相比其實沒有什麼變化，所以就有

\[h(p[j]*i)=h(i)\]

於是現在的答案變成了

\[Ans=\sum_{d=1}^N\frac{(\left \lfloor \frac{N}{d}\right \rfloor+1)\times \left \lfloor \frac{N}{d} \right \rfloor\times(\left \lfloor \frac{M}{d}\right \rfloor+1)\times \left \lfloor \frac{M}{d} \right \rfloor}{4}\times d\times h(d)\]

於是直接求\(d\times h(d)\)的字首和就好了

程式碼

#include<iostream>
#include<cstring>
#include<cstdio>
#define re register
#define maxn 10000005
#define LL long long
#define min(a,b) ((a)<(b)?(a):(b))
#define max(a,b) ((a)>(b)?(a):(b))
const LL mod=20101009;
int n,m;
int f[maxn],p[maxn>>2];
LL pre[maxn];
int main()
{
    scanf("%d%d",&n,&m);
    if(n>m) std::swap(n,m);
    f[1]=pre[1]=1;
    for(re int i=2;i<=n;i++)
    {
        if(!f[i]) p[++p[0]]=i,pre[i]=(1-i+mod);
        for(re int j=1;j<=p[0]&&p[j]*i<=n;j++)
        {
            f[p[j]*i]=1;
            if(i%p[j]==0) 
            {
                pre[i*p[j]]=pre[i];
                break;
            }
            pre[p[j]*i]=pre[p[j]]*pre[i];
        }
    }
    for(re int i=1;i<=n;i++) pre[i]=(i*pre[i]%mod+pre[i-1])%mod;
    LL ans=0;
    for(re LL l=1,r;l<=n;l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        LL N=n/l,M=m/l;
        ans=(ans+((N+1)*N/2%mod)*((M+1)*M/2%mod)%mod*(pre[r]-pre[l-1]+mod)%mod)%mod;
    }
    printf("%lld\n",ans);
    return 0;
}