设$f(x)=\sum_{x|d}p(d)$。

则$ans=\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^n\mu(i)\mu(j)\mu(k)f(lcm(i,j))f(lcm(i,k))f(lcm(j,k))$。

转化成图论模型，$i$到$j$有边的条件是$\mu(i)\neq0,\mu(j)\neq0,lcm(i,j)\leq n$。

枚举square-free的$\gcd$，再枚举square-free的$lcm$，然后枚举$\frac{lcm}{\gcd}$的因子$a$，那么可以得到一对满足条件的数对$(a\times\gcd,\frac{lcm}{a})$。

这样可以不重不漏地枚举出所有边，然后三元环计数即可。

时间复杂度$O(n\log n\sqrt{n\log n})$。

 

#include
     
      #include
      
       using namespace std;const int N=100005,M=1166760,E=760745;int n,i,j,k,x,y,z,P[N],f[N],vis[N],mu[N],p[N],tot,ans,A,B;int g[N],v[M],nxt[M],ed,d[N];inline void read(int&a){char c;while(!(((c=getchar())>='0')&&(c<='9')));a=c-'0';while(((c=getchar())>='0')&&(c<='9'))(a*=10)+=c-'0';}inline void add(int x,int y){v[++ed]=y;nxt[ed]=g[x];g[x]=ed;}namespace Triple{int g[N],v[E],w[E],nxt[E],ed,st[N],en[N],m,q[M],val[M],tmp[N];inline void add(int x,int y,int z){  if(d[x]>d[y])swap(x,y);  v[++ed]=y;w[ed]=z;nxt[ed]=g[x];g[x]=ed;}void solve(){  for(i=1;i<=n;i++)if(g[i]){    st[i]=m+1;    for(j=g[i];j;j=nxt[j])tmp[q[++m]=v[j]]=w[j];    en[i]=m;    sort(q+st[i],q+m+1);    for(j=st[i];j<=m;j++)val[j]=tmp[q[j]];  }  for(i=1;i<=n;i++)for(j=g[i];j;j=nxt[j]){    k=v[j],x=st[i],y=st[k];    while(x<=en[i]&&y<=en[k])if(q[x]==q[y]){      z=q[x];      A+=mu[i]*mu[k]*mu[z]*w[j]*val[x]*val[y];      x++,y++;    }else q[x]
       
        =y)continue;    d[x]++,d[y]++;    B+=(mu[x]*f[y]+mu[y]*f[x])*f[j]*f[j];  }  for(i=1;i<=n;i++)if(mu[i])for(j=i;j<=n;j+=i)if(mu[j]&&mu[j/i])for(k=g[j/i];k;k=nxt[k]){    x=i*v[k],y=j/v[k];    if(x>=y)continue;    Triple::add(x,y,f[j]);  }  Triple::solve();  return printf("%d",(ans+A*6+B*3)&((1<<30)-1)),0;}