Link
设(f(i)=sumlimits_{j=1}^isumlimits_{k=1}^ioperatorname{lcm}(gcd(i,j),gcd(i,k)))，那么(ans=sumlimits_{i=1}^nf(i))。
根据直觉(f(n))肯定是积性函数。
[egin{aligned}
f(n)&=sumlimits_{i=1}^nsumlimits_{j=1}^noperatorname{lcm}(gcd(n,i),gcd(n,j))\
&=sumlimits_{d|n}sumlimits_{g|n}operatorname{lcm}(d,g)sumlimits_{i=1}^n[gcd(n,i)=d]sumlimits_{j=1}^n[gcd(n,j)=g]\
&=sumlimits_{d|n}sumlimits_{g|n}operatorname{lcm}(d,g)varphi(frac nd)varphi(frac ng)
end{aligned}
]
先考虑怎么拆(operatorname{lcm})。
[egin{aligned}
operatorname{lcm}(n,m)&=frac{nm}{gcd(n,m)}\
&=sumlimits_{d|gcd(n,m)}frac ndfrac mdd[gcd(frac nd,frac md)=1]\
&=sumlimits_{d|gcd(n,m)}sumlimits_{g|gcd(frac nd,frac md)}mu(g)g^2frac n{dg}frac m{dg}d
end{aligned}
]
然后把这个代进去。
[egin{aligned}
f(n)&=sumlimits_{d|n}sumlimits_{g|n}operatorname{lcm}(d,g)varphi(frac nd)varphi(frac ng)\
&=sumlimits_{x|n}sumlimits_{y|n}sumlimits_{u|gcd(x,y)}sumlimits_{v|gcd(frac xu,frac yu)}mu(v)uv^2frac x{uv}frac y{uv}varphi(frac nx)varphi(frac ny)\
&=sumlimits_{u|n}usumlimits_{v|frac nu}mu(v)v^2(sumlimits_{w|frac n{uv}}wvarphi(frac n{uvw}))^2
end{aligned}
]
也就是说(f=id imes(mucdot id_2) imes(varphi imes id)^2)。（幂在点积意义下）
因此(f)是积性函数，但是(sum f)看上去就不太能杜教筛，因此考虑zzt筛。
[egin{aligned}
f(p^e)&=sumlimits_{i=1}^{p^e}sumlimits_{j=1}^{p^e}operatorname{lcm}(gcd(i,p^e),gcd(j,p^e))\
&=sumlimits_{i=0}^esumlimits_{j=0}^ep^{max(i,j)}varphi(p^{e-i})varphi(p^{e-j})\
&=sumlimits_{i=0}^ep^ivarphi(p^{e-i})(varphi(p^{e-i})+2sumlimits_{j=0}^{i-1}varphi(p^{e-j}))\
&=p^e(1+2sumlimits_{i=0}^{e-1}p^{e-i}-p^{e-i-1})+sumlimits_{i=0}^{e-1}(p-1)^2p^{e-1}(p^{e-i-1}+2sumlimits_{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
",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();
}