Description：

求 $\begin{align}\sum_{i=1}^n\sum_{j=1}^n\varphi(gcd(\varphi(i),\varphi(j)))\end{align}$

$1\le n \le 2\times 10^6$

Solution：

先枚举 $gcd$ ， $\begin{align}原式&=\sum_{d=1}^n\varphi(d)\sum_{i=1}^n\sum_{j=1}^n[gcd(\varphi(i),\varphi(j))=d]\end{align}$

然后似乎没什么办法了，然而可以枚举 $\varphi(i)$ 和 $\varphi(j)$ 。

$\begin{align}原式&=\sum_{d=1}^n\varphi(d)\sum_{i=1}^n\sum_{j=1}^ns[i]\times s[j]\times [gcd(i,j)=d]\end{align}$

$s[i]$ 表示 $[1,n]$ 中有多少个数的欧拉函数值为 $i$ ，后面就可以反演了。

$\begin{align}原式&=\sum_{d=1}^n\varphi(d)\sum_{i=1}^{\lfloor \frac n d\rfloor}\sum_{j=1}^{\lfloor \frac n d\rfloor}s[id]\times s[jd]\times [gcd(i,j)=1]\end{align}$

$\begin{align}原式&=\sum_{d=1}^n\varphi(d)\sum_{i=1}^{\lfloor \frac n d\rfloor}\sum_{j=1}^{\lfloor \frac n d\rfloor}s[id]\times s[jd]\times \sum_{k|i,k|j}\mu(k)\end{align}$

$\begin{align}原式&=\sum_{d=1}^n\varphi(d)\sum_{k=1}^{\lfloor\frac n d\rfloor}\mu(k)(\sum_{k|i,i=1}^{\lfloor \frac n d\rfloor}s[id])^2\end{align}$

$\begin{align}原式&=\sum_{d=1}^n\varphi(d)\sum_{k=1}^{\lfloor\frac n d\rfloor}\mu(k)(\sum_{i=1}^{\lfloor \frac n {dk}\rfloor}s[idk])^2\end{align}$

设 $\begin{align}f(k)=\sum_{i=1}^{\lfloor\frac n k\rfloor}s[ik]\end{align}$ ，则 $\begin{align}原式&=\sum_{d=1}^n\varphi(d)\sum_{k=1}^{\lfloor\frac n d\rfloor}\mu(k)f(dk)^2\end{align}$

$f$ 可以暴力预处理，因为由调和级数得 $\begin{align}\sum_{i=1}^n\frac 1 i=\ln n+\gamma\end{align}$

最后答案也可以暴力求。

枚举 $\varphi$ 是个好技巧。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cstring>

using namespace std;

int n;

#define MAXN 2000010

bool isprime[MAXN];

int prime[MAXN],tot = 0;

int mu[MAXN],phi[MAXN];

typedef long long ll;

int s[MAXN];

ll f[MAXN];

inline void work()

{

	scanf("%d",&n);

	tot = 0;

	for(register int i = 2;i <= n;++i)isprime[i] = true;

	mu[1] = 1;phi[1] = 1;

	for(register int i = 2;i <= n;++i)

	{

		if(isprime[i])

		{

			prime[++tot] = i;

			mu[i] = -1;phi[i] = i - 1;

		}

		for(register int j = 1;j <= tot && i * prime[j] <= n;++j)

		{

			register int k = i * prime[j];

			isprime[k] = false;

			if(i % prime[j] == 0)

			{

				mu[k] = 0;phi[k] = phi[i] * prime[j];

				break;

			}

			else

			{

				mu[k] = -mu[i];phi[k] = phi[i] * phi[prime[j]];

			}

		}

	}

	memset(s,0,sizeof(s));

	memset(f,0,sizeof(f));

	for(register int i = 1;i <= n;++i)++s[phi[i]];

	for(register int i = 1;i <= n;++i)

		for(register int j = i;j <= n;j += i)

			f[i] += s[j];

	register ll ans = 0;

	for(register int i = 1;i <= n;++i)

		for(register int j = 1;j <= n / i;++j)

			ans += phi[i] * mu[j] * f[i * j] * f[i * j];

	cout << ans << endl;

	return;

}

int main()

{

	int testcases = 0;

	scanf("%d",&testcases);

	while(testcases--)work(); 

	return 0;

}