原题链接

题目大意

给出两个整数$n,m(1\leq n,m \leq 10^7)$
$$\sum_{i=1}^n\sum_{j=1}^m lcm(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)}$$
枚举$gcd(i,j)$
$$\sum_{d=1}^{min(n,m)}\frac{1}{d}\sum_{i=1}^n\sum_{j=1}^m {ij}(gcd(i,j)=d)=\sum_{d=1}^{min(n,m)}\frac{1}{d}\sum_{i=1}^n\sum_{j=1}^m {ij}(gcd(\frac id,\frac jd)=1)$$
$$\sum_{d=1}^{min(n,m)}\frac{1}{d}\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor} {id jd}(gcd(i,j)=1)=\sum_{d=1}^{min(n,m)}d\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor} {ij}(gcd(i,j)=1)$$
考虑设$$sum(n,m)=\sum_{i=1}^n\sum_{j=1}^mij(gcd(i,j)=1)$$
根据莫比乌斯函数性质
$$sum(n,m)=\sum_{i=1}^n\sum_{j=1}^mij\sum_{d|gcd(i,j)}\mu(d)$$
枚举d可得
$$sum(n,m)=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}idjd=\sum_{d=1}^{min(n,m)}\mu(d)d^2\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}ij$$
不妨再设
$$G(n,m)=\sum_{i=1}^n\sum_{j=1}^m ij=\sum_{i=1}^n i\sum_{j=1}^m j=\frac{n(n+1)}{2}×\frac{m(m+1)}{2}$$
然后前缀和预处理$\mu(d)d^2$+数论分块即可

代码

#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=10000010;
const ll mod=20101009;
int prime[N],cnt;
bool st[N];
int mob[N];
int s[N];
int n,m;
void init(int n)
{
    mob[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!st[i]) prime[++cnt]=i,mob[i]=-1;

        for(int j=1;prime[j]<=n/i;j++)
        {
            st[prime[j]*i]=1;
            if(i%prime[j]==0)
            {
                mob[i*prime[j]]=0;
                break;
            }
            mob[i*prime[j]]=mob[i]*mob[prime[j]];
        }
    }
    for(int i=1;i<=n;i++) s[i]=(s[i-1]+1ll*mob[i]*i%mod*i%mod+mod)%mod;

}
int G(int n,int m)
{
    return (1ll*n*(n+1)/2%mod)*(1ll*m*(m+1)/2%mod)%mod;
}
int sum(int n,int m)
{
    int res=0;
    for(int l=1,r;l<=min(n,m);l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        res=(res+1ll*(s[r]-s[l-1])*G(n/l,m/l)%mod)%mod;
    }
    return (res+mod)%mod;
}
int solve(int n,int m)
{
    int res=0;
    for(int l=1,r;l<=min(n,m);l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        res=(res+1ll*(l+r)*(r-l+1)/2%mod*sum(n/l,m/l)%mod)%mod;
    }
    return (res+mod)%mod;
}
int main()
{
    int T=1;
    //cin>>T;
    while(T--)
    {
        cin>>n>>m;
        init(max(n,m));
        cout<<solve(n,m)<<'\n';
    }

    return 0;
}

要加油哦~