首先题目中给出的代码打错了，少了个等于号，应该是
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} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac{n}{d} \right \rfloor}[gcd(i,j)==d]ijd
\]
\]
然后需要一个打表找规律，发现\( \sum_{i=1}{n}\sum_{j=1}{i}[gcd(i,j)1]ij=\sum_{i=1}^{n}i\frac{i\phi(i)+[i1]}{2} \)
\[=\sum_{d=1}^{n}d((2)\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}i\frac{i\phi(i)+[i==1]}{2})-1)
\]
\]
\[=\sum_{d=1}^{n}d((\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}i^2\phi(i)+[i==1])-1)
\]
\]
\[=\sum_{d=1}^{n}d\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}i^2\phi(i)
\]
\]
然后就可以递归使用杜教筛了，关于用杜教筛求\( \sum_{i=1}{n}i2\phi(i) \)的前缀和，有如下推导：
设
\[g(n)=\sum_{i=1}^{n}i^2\sum_{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^2\phi(i)
$$那么把g展开：
\]
$$那么把g展开：
\]
g(n)=\sum_{i=2}{n}i2\sum_{d=1}^{i-1}[d|i]\phi(d)+s(n)
\[\]
s(n)=g(n)-\sum_{i=2}{n}i2\sum_{d=1}^{i-1}[d|i]\phi(d)
\[\]
=g(n)-\sum_{k=2}{n}k2\sum_{d=1}^{\left \lfloor \frac{n}{k} \right \rfloor}d\phi(d)
\[\]
=g(n)-\sum_{k=2}{n}k2*s(\left \lfloor \frac{n}{k} \right \rfloor)
\[\]
=\frac{n2(n+1)2}{4}-\sum_{k=2}{n}k2*s(\left \lfloor \frac{n}{k} \right \rfloor)
\[这就是标准的杜教筛递归子问题形式了，直接求解即可。
```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\n",(ans%mod+mod)%mod);
return 0;
}
```\]
```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\n",(ans%mod+mod)%mod);
return 0;
}
```\]