本题可能有更简单的做法,下面是我自己的做法:
可能需要的一些前置知识:



下面是推导:

for(int i=1;i<=n;i++){
    for(int j=1;j<=n;j++){
        ans[i]+=mu[j]*S(j,i/j)*S(j,i/j);
    }
}

细心观察之后就能发现,对于一个来说他能贡献给区间一个相同的值，那么做法就很显然了，预处理出和把所有答案算出来就能回答。

细节见代码:

#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 5e4 + 10;
#define fi first
#define se second
#define pb push_back
#define wzh(x) cerr<<#x<<'='<<x<<endl;
int a[N],cnt,p[N],mu[N];
LL phi[N];
void P(){
  phi[1]=mu[1]=1;
  for(int i=2;i<N;i++){
    if(!p[i])a[++cnt]=i,mu[i]=-1,phi[i]=i-1;
    for(int j=1;j<=cnt&&1ll*a[j]*i<N;j++){
      p[a[j]*i]=1;
      if(i%a[j]==0){
        mu[i*a[j]]=0;
        phi[i*a[j]]=a[j]*phi[i];
        break;
      }
      mu[a[j]*i]=-mu[i];
      phi[i*a[j]]=phi[i]*phi[a[j]];
    }
  }
}
vector<int>v[50003];
int L[N];
int main() {
  ios::sync_with_stdio(false);
  P();int cn=0;
  for(int i=1;i<=50000;i++){
    for(int j=1;j<=50000/i;j++){
      v[i].pb(mu[j*i]);
    }
    for(int j=1;j<v[i].size();j++){
      v[i][j]+=v[i][j-1];
    }
  }
  int t,n;
  for(int i=1;i<=50000;i++){
    for(int j=0,k=min(i-1,50000);j<=50000;){
      int no=k/i;
      if(no-1<(int)v[i].size()){
        assert(no-1<(int)v[i].size());
        int d=(no>=1)?(v[i][no-1]*v[i][no-1]):0;
        int now=mu[i]*d;
        assert(k+1<=50001);
        L[j]+=now;L[k+1]-=now;
      }
      j=k+1;
      k+=i;
      k=min(k,50000);
    }
  }
  for(int i=2;i<=50000;i++)L[i]+=L[i-1];
  for(cin>>t;t;t--){
    cin>>n;
    cout<<L[n]<<'\n';
  }
  return 0;
}