首先由这样一个结论：
[d(ij)=sum_{p|i}sum_{q|j}[gcd(p,q)==1]
]
然后推反演公式：
[sum_{i=1}^{n}sum_{j=1}^{n}sum_{p|i}sum_{q|j}[gcd(p,q)==1]
]
[sum_{p=1}^{n}sum_{q=1}^{n}[gcd(p,q)==1]left lfloor frac{n}{p}
ight
floorleft lfloor frac{n}{q}
ight
floor
]
[sum_{p=1}^{n}sum_{q=1}^{n}sum_{k|p,k|q}mu(k)left lfloor frac{n}{p}
ight
floorleft lfloor frac{n}{q}
ight
floor
]
[sum_{k=1}^{n}mu(k)(sum_{k|p}left lfloor frac{n}{p}
ight
floor)^2
]
[sum_{k=1}^{n}mu(k)(sum_{p=1}^{left lfloor frac{n}{k}
ight
floor}left lfloor frac{n}{pk}
ight
floor)^2
]
然后对于这个递归子问题形式就可以用杜教筛求解了
#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
const long long N=2000005,m=2000000,P=100005,mod=1e9+7;
long long n,mb[N],q[N],tot,ans,p[P],t;
bool v[N];
long long F(long long n)
{
long long sum=0;
for(long long i=1,la;i<=n;i=la+1)
{
la=n/(n/i);
sum=(sum+(la-(i-1))*(n/i)%mod)%mod;
}
return sum*sum%mod;
}
long long getp(long long x,long long n)
{
return (x<=m)?mb[x]:p[n/x];
}
void slv(long long x,long long n)
{
if(x<=m)
return;
long long i,j=1,t=n/x;
if(v[t])
return;
v[t]=1;
p[t]=1;
while(j<x)
{
i=j+1;
j=x/(x/i);
slv(x/i,n);
p[t]=(p[t]-getp(x/i,n)*(j-i+1)+mod)%mod;
}
}
long long wk(long long n)
{
if(n<=m)
return mb[n];
memset(v,0,sizeof(v));
slv(n,n);
return p[1];
}
int main()
{
mb[1]=1;
for(long long i=2;i<=m;i++)
{
if(!v[i])
{
q[++tot]=i;
mb[i]=-1;
}
for(long long j=1;j<=tot&&i*q[j]<=m;j++)
{
long long k=i*q[j];
v[k]=1;
if(i%q[j]==0)
{
mb[k]=0;
break;
}
mb[k]=-mb[i];
}
}
for(long long i=1;i<=m;i++)
mb[i]+=mb[i-1];
scanf("%lld",&n);
for(long long i=1,la;i<=n;i=la+1)
{
la=n/(n/i);//cout<<la<<" "<<i-1<<" "<<wk(la)<<" "<<wk(i-1)<<endl;
ans=(ans+(wk(la)-wk(i-1)+mod)%mod*F(n/i)%mod)%mod;
}
printf("%lld
",ans);
return 0;
}