首先由這樣一個結論：
\[ 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} \right \rfloor\left \lfloor \frac{n}{q} \right \rfloor \]

#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\n",ans);
    return 0;
}

bzoj 4176: Lucas的數論【莫比烏斯反演+杜教篩】