题面
题目分析
（默认\(n<m\)）
题目要求\(\sum\limits_{i=1}^n\sum\limits_{j=1}^mlcm(i,j)\)。
由\(lcm(i,j)=\frac{i\cdot j}{gcd(i,j)}\)
得:
\[ \begin{split} ans & =\sum\limits_{i=1}^n\sum\limits_{j=1}^m\frac{i\cdot j}{gcd(i,j)} \\ & =\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^m\frac{i\cdot j}{d}[gcd(i,j)==d]\\ & = \sum\limits_{d=1}^nd\cdot \sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}i\cdot j[gcd(i,j)==1] \end{split} \]
如果只看最后一部分，\(\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}i\cdot j[gcd(i,j)==1]\)，可以很自然想到莫比乌斯反演。
（以下\(n,m,gcd\)为\(\lfloor\frac n d\rfloor,\lfloor\frac m d \rfloor,\lfloor\frac {gcd}d\rfloor\)）
我们用\(g(i)\)表示\(gcd(i,j)==k\cdot i,k\in Z\)的贡献，\(f(i)\)表示\(gcd(i,j)==i\)的贡献。
于是有\(g(x)=\sum\limits_{x|d}^nf(d) \Rightarrow f(x)=\sum\limits_{x|d}^n\mu(\frac d x)\cdot g(d)\)。
只要可以快速求出\(g(d)\)便可得到答案。
对于\(g(x)\)：
\[ \begin{split} g(x)&=\sum\limits_{i=1}^n\sum\limits_{j=1}^mi\cdot j[x|gcd(i,j)]\\ &=x\cdot x\cdot \sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{x}\rfloor}i\cdot j[1|gcd(i,j)]\\ &=x\cdot x\cdot \sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}i\sum\limits_{j=1}^{\lfloor\frac{m}{x}\rfloor}j\\ &=x\cdot x\cdot \frac{(1+\lfloor\frac{n}{x}\rfloor)\cdot \lfloor\frac{n}{x}\rfloor}{2}\cdot \frac{(1+\lfloor\frac{m}{x}\rfloor)\cdot \lfloor\frac{m}{x}\rfloor}{2} \end{split} \]最终
\[ \begin{split} ans &=\sum\limits_{d=1}^n d\cdot f(1)\\ &=\sum\limits_{d=1}^n d\cdot \sum\limits_{i=1}^{\lfloor\frac n d\rfloor}\mu(i)\cdot g(i)\\ \end{split} \]你会发现，现在的时间复杂度还是有问题，这时候就需要整除分块求解。
P.S
加强版：
代码实现
#include#include #include #include #include #include #include #define MAXN 0x7ffffffftypedef 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;}