解法：莫比乌斯反演

更新了更优化的算法，在最下面…

设$f(d) = [gcd = d]$，表示gcd为d的数的对数的个数

设$g(p) = \sum_{p|d}f(d)$，表示gcd为p的倍数的数对的个数，其中$p$代表的质数

易知$g(p) = \lfloor n / p \rfloor \lfloor n / p \rfloor $ ，因为如果两个数都是p的倍数，那他们的gcd也一定是p的倍数。小于n的p的倍数的个数就是$\lfloor n / p \rfloor$，由于需要的是两个组成一个数对，所以乘法原理相乘一下。

则$g(p) = \sum_{p|d}f(d) = \lfloor n / p \rfloor \lfloor n / p \rfloor$里的g是好求的，而我们需要的是f函数和的值，可以用莫比乌斯反演，将计数关系交换。

反演以后$f(p) = \sum_{p|d} \mu(d/p) g(d) $

将函数g展开即$f(p) = \sum_{p|d} \mu(d/p) \lfloor n / d \rfloor \lfloor n / d \rfloor $

其中mu为莫比乌斯函数，这一项可以在线筛的时候求和，后面下取整除法只有$2\sqrt{n}$种取法，所以可以分块计算。

计算一个质数的f函数复杂度$O(\sqrt{N})$，题目给出的N为1e7，根据素数定理，小于N的质数个数约为$N / logN$个，也就是大概3e5个。一个一个算复杂度大概3e8可以接受

总复杂度$O(N / logN * \sqrt{n} )$

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 1e7 + 233;
int primes[maxn], mu[maxn], sum[maxn], cnt;
bool st[maxn];
void get_primes(int n)
{
    mu[1] = 1; sum[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!st[i]) primes[cnt++] = i, mu[i] = -1;
        sum[i] = sum[i - 1] + mu[i];
        for(int j = 0; j < cnt &&  i * primes[j] <= n; j++)
        {
            st[primes[j] * i] = 1;
            if(i % primes[j] == 0)
            {
                mu[primes[j] * i] = 0;
                break;
            }
            else mu[primes[j] * i] = -mu[i];
        }
    }
}
int main()
{
    get_primes(10000000);
    int n; cin >> n;
    ll ans = 0;
    for(int j = 0; j < cnt && primes[j] <= n; j++)
    {
        int a = n / primes[j], c = 0;
        for(int i = 1; i <= a; i = c + 1)
        {
            c = n / (n / i);
            ll b = i * primes[j];
            ll t = (n / b) * (n / b);
            ans += (ll)(sum[c] - sum[i - 1]) * t;
        }
    }
    cout<<ans;
}

更新

注意到$f(p) = \sum_{p|d} \mu(d/p) g(d) $，和答案 $ ans = \sum _ {p} f(p)$ 的关系
即:
$ \sum _ {p} \sum_{p|d} \mu(d/p) \lfloor n / d \rfloor \lfloor n / d \rfloor $

$ \sum _ {p} \sum_{d = 1} \mu(d) \lfloor n / pd \rfloor \lfloor n / pd \rfloor $

考虑把除法分块的部分拿到外面来，可以减少一层循环。
$pd$的变量与内层外层都有关，要拿到外面去要去掉内层d变量的影响，于是设$T = pd$，由于T是质数的倍数，所以可以取到除1以外的所有数字。
同时内层变量被取代，$d = T / p$

$ \sum _ {T = 1} \lfloor n / T \rfloor \lfloor n / T \rfloor \sum _ {p|T}\mu{(T/p)}$

有了这个式子，设$h(T) = \sum _ {p|T}\mu{(T/p)}$ ，想办法预处理它。可以对于每个质数，将它的倍数的函数值加上倍数

前面提到1~n的质数个数约$n / \log(n)$个，对于每个质数，处理到n需要的复杂度约$log(n)$，所以预处理f函数的均摊复杂度接近$O(n)$

然后再分块处理出答案，这部分只有$\sqrt{n}$.

总复杂度$O(n + \sqrt{n})$。

本以为可以AC的时候，惊人的MLE出现了.....，于是考虑到mu只有-1，0，1三种值，我把它改成了short....

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 1e7 + 233;
int primes[maxn], f[maxn], sum[maxn], cnt;
short mu[maxn];
bool st[maxn];
void get_primes(int n)
{
    mu[1] = 1; 
    for(int i = 2; i <= n; i++)
    {
        if(!st[i]) primes[cnt++] = i, mu[i] = -1;
        for(int j = 0; j < cnt &&  i * primes[j] <= n; j++)
        {
            st[primes[j] * i] = 1;
            if(i % primes[j] == 0)
            {
                mu[primes[j] * i] = 0;
                break;
            }
            else mu[primes[j] * i] = -mu[i];
        }
    }
}
int main()
{
    get_primes(10000000);
    int n; cin >> n;
    for(int i = 0; i < cnt && primes[i] <= n; i++)
    {
        for(int j = 1; j <= n / primes[i]; j++)
        {
            f[primes[i] * j] += mu[j];
        }
    }
    for(int i = 1; i <= n; i++) sum[i] = sum[i - 1] + f[i];
    ll ans = 0;

    int c = 0;
    for(int i = 2; i <= n; i = c + 1)
    {
        c = n / (n / i);
        ll t = (ll)(n / i) * (n / i);
        ans += (ll)(sum[c] - sum[i - 1]) * t;
    }

    cout<<ans;
}