题面
题目分析
（默认(n<m)）
题目要求(sumlimits_{i=1}^nsumlimits_{j=1}^mlcm(i,j))。
由(lcm(i,j)=frac{icdot j}{gcd(i,j)})
得:
[egin{split}
ans & =sumlimits_{i=1}^nsumlimits_{j=1}^mfrac{icdot j}{gcd(i,j)} \
& =sumlimits_{d=1}^nsumlimits_{i=1}^nsumlimits_{j=1}^mfrac{icdot j}{d}[gcd(i,j)==d]\
& = sumlimits_{d=1}^ndcdot sumlimits_{i=1}^{lfloorfrac{n}{d}
floor}sumlimits_{j=1}^{lfloorfrac{m}{d}
floor}icdot j[gcd(i,j)==1]
end{split}
]
如果只看最后一部分，(sumlimits_{i=1}^{lfloorfrac{n}{d} floor}sumlimits_{j=1}^{lfloorfrac{m}{d} floor}icdot j[gcd(i,j)==1])，可以很自然想到莫比乌斯反演。
（以下(n,m,gcd)为(lfloorfrac n d floor,lfloorfrac m d floor,lfloorfrac {gcd}d floor)）
我们用(g(i))表示(gcd(i,j)==kcdot i,kin Z)的贡献，(f(i))表示(gcd(i,j)==i)的贡献。
于是有(g(x)=sumlimits_{x|d}^nf(d) Rightarrow f(x)=sumlimits_{x|d}^nmu(frac d x)cdot g(d))。
只要可以快速求出(g(d))便可得到答案。
对于(g(x))：
[egin{split}
g(x)&=sumlimits_{i=1}^nsumlimits_{j=1}^micdot j[x|gcd(i,j)]\
&=xcdot xcdot sumlimits_{i=1}^{lfloorfrac{n}{x}
floor}sumlimits_{j=1}^{lfloorfrac{m}{x}
floor}icdot j[1|gcd(i,j)]\
&=xcdot xcdot sumlimits_{i=1}^{lfloorfrac{n}{x}
floor}isumlimits_{j=1}^{lfloorfrac{m}{x}
floor}j\
&=xcdot xcdot frac{(1+lfloorfrac{n}{x}
floor)cdot lfloorfrac{n}{x}
floor}{2}cdot frac{(1+lfloorfrac{m}{x}
floor)cdot lfloorfrac{m}{x}
floor}{2}
end{split}
]
最终
[egin{split}
ans &=sumlimits_{d=1}^n dcdot f(1)\
&=sumlimits_{d=1}^n dcdot sumlimits_{i=1}^{lfloorfrac n d
floor}mu(i)cdot g(i)\
end{split}
]
你会发现，现在的时间复杂度还是有问题，这时候就需要整除分块求解。
P.S
加强版：【BZOJ2693】jzptab
代码实现
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<iomanip>
#include<cstdlib>
#define MAXN 0x7fffffff
typedef long long LL;
const int N=1e7+5,mod=20101009;
using namespace std;
inline int Getint(){register int x=0,f=1;register char ch=getchar();while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}return x*f;}
int mu[N],prime[N],g[N];
bool vis[N];
int t(int x){return 1ll*x*(x+1)/2%mod;}
int f(int n,int m){
if(n>m)swap(n,m);
LL ans=0;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+1ll*(g[r]-g[l-1])*t(n/l)%mod*t(m/l)%mod)%mod;
}
return ans;
}
int main(){
mu[1]=g[1]=1;
for(int i=2;i<=1e7;i++){
if(!vis[i])prime[++prime[0]]=i,mu[i]=-1;
for(int j=1;j<=prime[0]&&1ll*i*prime[j]<=1e7;j++){
vis[i*prime[j]]=1;
if(i%prime[j]==0)break;
mu[i*prime[j]]=-mu[i];
}
g[i]=(g[i-1]+1ll*i*i*mu[i]%mod)%mod;
}
int n=Getint(),m=Getint();
if(n>m)swap(n,m);
LL ans=0;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+1ll*(l+r)*(r-l+1)/2%mod*f(n/l,m/l)%mod)%mod;
}
cout<<(ans+mod)%mod;
return 0;
}