首先由这样一个式子：( d(ij)=sum_{p|i}sum_{q|j}[gcd(p,q)==1]frac{pj}{q} )大概感性证明一下吧我不会证
然后开始推：
[sum_{i=1}^{n}sum_{j=1}^{n}sum_{p|i}sum_{q|j}[gcd(p,q)==1]frac{pj}{q}
]
[sum_{p=1}^{n}sum_{q=1}^{n}[gcd(p,q)==1]sum_{p|i}sum_{q|j}frac{pj}{q}
]
[sum_{p=1}^{n}psum_{q=1}^{n}[gcd(p,q)==1]left lfloor frac{n}{p}
ight
floorsum_{j=1}^{left lfloor frac{n}{q}
ight
floor}j
]
方便起见设( f(n)=sum_{i=1}^{n}i )
[sum_{p=1}^{n}psum_{q=1}^{n}sum_{k|p,k|q}mu(k)left lfloor frac{n}{p}
ight
floor f(left lfloor frac{n}{q}
ight
floor)
]
[sum_{k=1}^{n}mu(k)sum_{k|p}pleft lfloor frac{n}{p}
ight
floorsum_{k|q}f(left lfloor frac{n}{q}
ight
floor)
]
[sum_{k=1}^{n}mu(k)sum_{i=1}^{left lfloor frac{n}{k}
ight
floor}ikleft lfloor frac{n}{ik}
ight
floorsum_{j=1}^{left lfloor frac{n}{k}
ight
floor}f(left lfloor frac{n}{jk}
ight
floor)
]
[sum_{k=1}^{n}mu(k)ksum_{i=1}^{left lfloor frac{n}{k}
ight
floor}ileft lfloor frac{n}{ik}
ight
floorsum_{j=1}^{left lfloor frac{n}{k}
ight
floor}f(left lfloor frac{n}{jk}
ight
floor)
]
这个样子显然可以用杜教筛了，但是注意到后面有两个求和式，可能会增大常数（但是也不会T啦），所以考虑这两个求和式的关系：
[sum_{i=1}^{n}f(left lfloor frac{n}{i}
ight
floor)
]
[=sum_{i=1}^{n}sum_{j=1}^{left lfloor frac{n}{i}
ight
floor}j
]
[=sum_{i=1}^{n}sum_{j=1}^{left lfloor frac{n}{i}
ight
floor}j
]
[=sum_{j=1}^{n}jleft lfloor frac{n}{j}
ight
floor
]
所以这两个式子是一样的！于是就变成了：
[sum_{k=1}^{n}mu(k)k(sum_{j=1}^{left lfloor frac{n}{k}
ight
floor}f(left lfloor frac{n}{jk}
ight
floor))^2
]
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int N=1000005,inv2=500000004,mod=1e9+7;
int m,mb[N],q[N],tot;
long long n,s[N],ans,ha[N];
bool v[N];
long long slv(long long n)
{
return n*(n+1)%mod*inv2%mod;
}
long long wk(long long x)
{
if(x<=m)
return s[x];//cout<<x<<endl;
if(ha[n/x])
return ha[n/x];
long long re=1ll;
for(int i=2,la;i<=x;i=la+1)
{
la=x/(x/i);
re=(re-(slv(la)-slv(i-1))*wk(x/i)%mod)%mod;
}
return ha[n/x]=re;
}
long long clc(long long n)
{
long long re=0ll;
for(int i=1,la;i<=n;i=la+1)
{
la=n/(n/i);
re=(re+(la-i+1)*slv(n/i)%mod)%mod;
}
return re;
}
int main()
{
scanf("%lld",&n);
m=(int)ceil(pow((int)n,2.0/3));
mb[1]=1;
for(int i=2;i<=m;i++)
{
if(!v[i])
{
q[++tot]=i;
mb[i]=-1;
}
for(int j=1;j<=tot&&q[j]*i<=m;j++)
{
int k=i*q[j];
v[k]=1;
if(i%q[j]==0)
{
mb[k]=0;
break;
}
mb[k]=-mb[i];
}
}
for(int i=1;i<=m;i++)
s[i]=(s[i-1]+i*mb[i])%mod;
//cout<<wk(102)<<" "<<wk(101)<<endl;
for(int i=1,la;i<=n;i=la+1)
{
la=n/(n/i);
long long ml=clc(n/i);//if(i!=1)cout<<i-1<<" "<<n/(i-1)<<endl<<la<<" "<<n/la<<endl;
ans=(ans+(wk(la)-wk(i-1))*ml%mod*ml%mod)%mod;
}
printf("%lld",(ans%mod+mod)%mod);
return 0;
}