Link
设\(f(i)=\sum\limits_{j=1}^i\sum\limits_{k=1}^i\operatorname{lcm}(\gcd(i,j),\gcd(i,k))\)，那么\(ans=\sum\limits_{i=1}^nf(i)\)。
根据直觉\(f(n)\)肯定是积性函数。

\[\begin{aligned} f(n)&=\sum\limits_{i=1}^n\sum\limits_{j=1}^n\operatorname{lcm}(\gcd(n,i),\gcd(n,j))\\ &=\sum\limits_{d|n}\sum\limits_{g|n}\operatorname{lcm}(d,g)\sum\limits_{i=1}^n[\gcd(n,i)=d]\sum\limits_{j=1}^n[\gcd(n,j)=g]\\ &=\sum\limits_{d|n}\sum\limits_{g|n}\operatorname{lcm}(d,g)\varphi(\frac nd)\varphi(\frac ng) \end{aligned} \]

先考虑怎么拆\(\operatorname{lcm}\)。

\[\begin{aligned} \operatorname{lcm}(n,m)&=\frac{nm}{\gcd(n,m)}\\ &=\sum\limits_{d|\gcd(n,m)}\frac nd\frac mdd[\gcd(\frac nd,\frac md)=1]\\ &=\sum\limits_{d|\gcd(n,m)}\sum\limits_{g|\gcd(\frac nd,\frac md)}\mu(g)g^2\frac n{dg}\frac m{dg}d \end{aligned} \]

然后把这个代进去。

\[\begin{aligned} f(n)&=\sum\limits_{d|n}\sum\limits_{g|n}\operatorname{lcm}(d,g)\varphi(\frac nd)\varphi(\frac ng)\\ &=\sum\limits_{x|n}\sum\limits_{y|n}\sum\limits_{u|\gcd(x,y)}\sum\limits_{v|\gcd(\frac xu,\frac yu)}\mu(v)uv^2\frac x{uv}\frac y{uv}\varphi(\frac nx)\varphi(\frac ny)\\ &=\sum\limits_{u|n}u\sum\limits_{v|\frac nu}\mu(v)v^2(\sum\limits_{w|\frac n{uv}}w\varphi(\frac n{uvw}))^2 \end{aligned} \]

也就是说\(f=id\times(\mu\cdot id_2)\times(\varphi\times id)^2\)。（幂在点积意义下）
因此\(f\)是积性函数，但是\(\sum f\)看上去就不太能杜教筛，因此考虑zzt筛。

\[\begin{aligned} f(p^e)&=\sum\limits_{i=1}^{p^e}\sum\limits_{j=1}^{p^e}\operatorname{lcm}(\gcd(i,p^e),\gcd(j,p^e))\\ &=\sum\limits_{i=0}^e\sum\limits_{j=0}^ep^{\max(i,j)}\varphi(p^{e-i})\varphi(p^{e-j})\\ &=\sum\limits_{i=0}^ep^i\varphi(p^{e-i})(\varphi(p^{e-i})+2\sum\limits_{j=0}^{i-1}\varphi(p^{e-j}))\\ &=p^e(1+2\sum\limits_{i=0}^{e-1}p^{e-i}-p^{e-i-1})+\sum\limits_{i=0}^{e-1}(p-1)^2p^{e-1}(p^{e-i-1}+2\sum\limits_{j=0}^{i-1}p^{e-j-1})\\ &=(2e+1)(p^{2e}-p^{2e-1})+p^{e-1} \end{aligned} \]

\(f(p)=3p^2-3p+1\)。

#include<cstdio>
using u32=unsigned int;
using u64=unsigned long long;
const int N=100007;
int n,cnt,tot,lim,pr[N],is[N],num[N];u32 h[N];
int read(){int x;scanf("%d",&x);return x;}
int getid(int x){return x<=lim? x:cnt+1-n/x;}
u64 S(u64 p,int deg){return !deg? p:deg==1? p*(p+1)/2:p%6==1? (p+1)*(p+p+1)/6*p:p%6==4? p*(p+p+1)/6*(p+1):p*(p+1)/6*(p+p+1);}
u32 g(int n,int m)
{
    if(n<pr[m]) return 0;
    u32 ans=h[getid(n)]-h[pr[m]];u64 p,q;
    for(int i=m+1,e;i<=tot&&pr[i]*pr[i]<=n;++i) for(p=pr[i],e=1,q=p*(p-1);p<=(u64)n;++e,p*=pr[i],q*=(u64)pr[i]*pr[i]) ans+=((e+e+1)*q+p/pr[i])*(g(n/p,i)+(e>1));
    return ans;
}
void init()
{
    static u32 t[3][N];cnt=0;
    for(lim=1;(lim+1)*(lim+1)<=n;++lim);
    for(int l=1,r;l<=n;l=r+1) r=n/(n/l),num[++cnt]=r;
    for(int i=1;i<=cnt;++i) for(int k=0;k<3;++k) t[k][i]=S(num[i],k)-1;
    for(int i=1;i<=tot&&pr[i]<=lim;++i) for(int j=cnt;j&&pr[i]*pr[i]<=num[j];--j) for(int k=0;k<3;++k) t[k][j]-=(k>=1? pr[i]:1)*(k>=2? pr[i]:1)*(t[k][getid(num[j]/pr[i])]-t[k][pr[i-1]]);
    for(int i=1;i<=cnt;++i) h[i]=3*(t[2][i]-t[1][i])+t[0][i];
}
void solve()
{
    n=read();
    init();
    printf("%u\n",g(n,0)+1);
}
int main()
{
    for(int i=2,j;i<=40000;++i) if(!is[i]) for(pr[++tot]=i,j=i*2;j<=40000;j+=i) is[j]=1;
    for(int T=read();T;--T) solve();
}