题目链接 数论模板 求$\sum_{i=1}^n\sum_{j=i+1}^n\gcd(i,j)$。

解法一：蓝书解法，计算贡献

设$f(n)=\sum_{i=1}^{n-1}\gcd(i,n)$，则所求答案为f(n)的前缀和。
注意到$\gcd(i,n)=k$的值一定是n的约数，按照这个约数进行分类。设满足$\gcd(x,n)=k$的约束有g(n,i)个，则$f(n)=sum{i\times g(n,i)|d\mod i==0}$。
依次计算f(n)需要枚举每个n的约数i，速度较慢，但是对每个i枚举它的倍数ki来更新f(ki)的值，此时速度和普通筛法同阶。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
struct EulerSieve
{
	vector<int> p,m,phi;
	EulerSieve(int N):m(N,0),phi(N,0)
	{
		phi[1]=1;
		for(long long i=2,k; i<N; ++i)
		{
			if(!m[i])p.push_back(m[i]=i),phi[i]=i-1;
			for(int j=0; j<p.size()&&(k=i*p[j])<N; ++j)
			{
				phi[k]=phi[i]*p[j];
				if((m[k]=p[j])==m[i])break;
				phi[k]-=phi[i];
			}
		}
	}
};
struct GCDExtremeII:EulerSieve
{
	vector<ll> f;
	GCDExtremeII(int N):EulerSieve(N),f(N)
	{
		for(int i=1; i<N; ++i)
		{
			for(int k=2; k*i<N; ++k)
				f[k*i]+=i*phi[k];
			f[i]+=f[i-1];
		}
	}
} g(4e6+9);
int main()
{
	for(int n; ~scanf("%d",&n)&&n;)printf("%lld\n",g.f[n]);
}

解法二：分块

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
struct EulerSieve
{
	vector<int> p,m,phi;
	vector<ll> sum;
	EulerSieve(int N):m(N,0),phi(N,0),sum(N,0)
	{
		phi[1]=1;
		for(long long i=2,k; i<N; ++i)
		{
			if(!m[i])p.push_back(m[i]=i),phi[i]=i-1;
			for(int j=0; j<p.size()&&(k=i*p[j])<N; ++j)
			{
				phi[k]=phi[i]*p[j];
				if((m[k]=p[j])==m[i])break;
				phi[k]-=phi[i];
			}
		}
		for(int i=0; i<N; ++i)sum[i]=sum[i-1]+phi[i];
	}
} e(4e6+9);
int main()
{
	for(ll n,sum,tmp; ~scanf("%lld",&n)&&n;)
	{
		for(int i=1,j=sum=0; i<=n; i=j+1)
			tmp=n/i,sum+=tmp*tmp*(e.sum[j=n/tmp]-e.sum[i-1]);
		printf("%lld\n",(sum-(n+1)*n/2)/2);
	}
}

解法三：莫比乌斯反演

$ans=\sum_{d=1}^nd\sum_{i=1}^{n/d}\mu(i)[\frac{n}{id}]^2=\sum_{T=1}^n[\frac{n}{T}]^2\sum_{d|T}d\mu(\frac{T}{d})$ 数论分块可$O(\sqrt n)$查询，总时间复杂度$O(n\sqrt n)$。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
struct EulerSieve
{
	vector<int> p,m,mu,sum;//素数序列，最小素因子，欧拉函数，莫比乌斯函数
	EulerSieve(int N):m(N,0),mu(N,0),sum(N,0)
	{
		mu[1]=1;//m[1]=0
		for(long long i=2,k; i<N; ++i)//防i*p[j]爆int
		{
			if(!m[i])p.push_back(m[i]=i),mu[i]=-1;//i是素数
			for(int j=0; j<p.size()&&(k=i*p[j])<N; ++j)
			{
				if((m[k]=p[j])==m[i])
				{
					mu[k]=0;
					break;
				}
				mu[k]=-mu[i];
			}
		}
		for(int i=1; i<N; ++i)sum[i]=sum[i-1]+mu[i];
	}
} e(4e6+9);
int main()
{
	for(ll n,sum,tmp; ~scanf("%lld",&n)&&n;)
	{
		sum=0;
		for(int i=1; i<=n; ++i)
			for(ll l=1,r,tmp,m=n/i; l<=m; l=r+1)
				tmp=m/l,sum+=i*tmp*tmp*(e.sum[r=m/tmp]-e.sum[l-1]);
		printf("%lld\n",(sum-(n+1)*n/2)/2);
	}
}