首先题目中给出的代码打错了，少了个等于号，应该是
G=0;
for(i=1;i<=N;i++)
for(j=1;j<=N;j++)
{
G = (G + lcm(i,j)) % 1000000007;
}
然后就是大力推公式：
[sum_{i=1}^{n}sum_{j=1}^{n}lcm(i,j)
]
[=sum_{i=1}^{n}sum_{j=1}^{n}frac{ij}{gcd(i,j)}
]
[=sum_{d=1}^{n}sum_{i=1}^{n}sum_{j=1}^{n}[gcd(i,j)==d]frac{ij}{d}
]
[=sum_{d=1}^{n}sum_{i=1}^{left lfloor frac{n}{d}
ight
floor}sum_{j=1}^{left lfloor frac{n}{d}
ight
floor}[gcd(i,j)==d]ijd
]
然后需要一个打表找规律，发现( sum_{i=1}{n}sum_{j=1}{i}[gcd(i,j)1]ij=sum_{i=1}^{n}ifrac{iphi(i)+[i1]}{2} )
[=sum_{d=1}^{n}d((2)sum_{i=1}^{left lfloor frac{n}{d}
ight
floor}ifrac{iphi(i)+[i==1]}{2})-1)
]
[=sum_{d=1}^{n}d((sum_{i=1}^{left lfloor frac{n}{d}
ight
floor}i^2phi(i)+[i==1])-1)
]
[=sum_{d=1}^{n}dsum_{i=1}^{left lfloor frac{n}{d}
ight
floor}i^2phi(i)
]
然后就可以递归使用杜教筛了，关于用杜教筛求( sum_{i=1}{n}i2phi(i) )的前缀和，有如下推导：
设
[g(n)=sum_{i=1}^{n}i^2sum_{d=1}^{i}[d|i]phi(d)=sum_{i=1}^{n}i^3=frac{n^2(n+1)^2}{4}
]
[s(n)=sum_{i=1}^{n}i^2phi(i)
$$那么把g展开：
]
g(n)=sum_{i=2}{n}i2sum_{d=1}^{i-1}[d|i]phi(d)+s(n)
[]
s(n)=g(n)-sum_{i=2}{n}i2sum_{d=1}^{i-1}[d|i]phi(d)
[]
=g(n)-sum_{k=2}{n}k2sum_{d=1}^{left lfloor frac{n}{k} ight floor}dphi(d)
[]
=g(n)-sum_{k=2}{n}k2*s(left lfloor frac{n}{k} ight floor)
[]
=frac{n2(n+1)2}{4}-sum_{k=2}{n}k2*s(left lfloor frac{n}{k} ight floor)
[这就是标准的杜教筛递归子问题形式了，直接求解即可。
```cpp
#include<iostream>
#include<cstdio>
using namespace std;
const long long N=1000005,m=1000000,inv2=500000004,inv4=250000002,inv6=166666668,mod=1e9+7;
long long n,phi[N],q[N],tot,ans,ha[N];
bool v[N];
long long wk1(long long x)
{
if(x>=mod)
x-=mod;
return x%mod*(x+1)%mod*inv2%mod;
}
long long wk2(long long x)
{
if(x>=mod)
x-=mod;
return x%mod*(x+1)%mod*(x%mod*2+1)%mod*inv6%mod;
}
long long wk3(long long x)
{
if(x>=mod)
x-=mod;
return x%mod*x%mod*(x+1)%mod*(x+1)%mod*inv4%mod;
}
long long slv(long long x)
{
if(x<=m)
return phi[x];
if(ha[n/x])
return ha[n/x];
long long re=wk3(x);
for(long long i=2,la;i<=x;i=la+1)
{
la=x/(x/i);
re=(re-(wk2(la)-wk2(i-1))%mod*slv(x/i)%mod)%mod;
}
return ha[n/x]=re;
}
int main()
{
phi[1]=1;
for(int i=2;i<=m;i++)
{
if(!v[i])
{
q[++tot]=i;
phi[i]=i-1;
}
for(int j=1;j<=tot&&i%mod*q[j]<=m;j++)
{
int k=i%mod*q[j];
v[k]=1;
if(i%q[j]==0)
{
phi[k]=phi[i]%mod*q[j];
break;
}
phi[k]=phi[i]%mod*(q[j]-1);
}
}
for(int i=1;i<=m;i++)
phi[i]=(phi[i]%mod*i%mod*i%mod+phi[i-1])%mod;
scanf("%lld",&n);
for(long long i=1,la;i<=n;i=la+1)
{
la=n/(n/i);
ans=(ans+(wk1(la)-wk1(i-1))%mod*slv(n/i)%mod)%mod;
}
printf("%lld
",(ans%mod+mod)%mod);
return 0;
}
```]