vector<int> prime; /// prime list
vector<int> lpf;   /// Lowest prime factor, lpf[x] is smallest prime divisor of x
void sieve(int lim = LIM) /// O(n)
{
    prime.assign(1, 2);
    lpf.assign(lim + 1, 2);

    lpf[1] = 1; /// For easier calculation but can cause inf loops
    for (int i = 3; i <= lim; i += 2) {
        if (lpf[i] == 2) prime.push_back(lpf[i] = i);
        for (int j = 0; j < sz(prime) && prime[j] <= lpf[i] && prime[j] * i <= lim; ++j)
            lpf[prime[j] * i] = prime[j];
    }
}

/// mask(x) is smallest positive number that mask(x) * x is a perfect square
int getMask(int x) /// O(log n)
{
    int mask = 1;
    while (x > 1) {
        int p = lpf[x], f = 0;
        do x /= p, f++; while (p == lpf[x]);
        if (f & 1) mask *= p; /// if current power is odd, we mutiple mask with current prime
    }
    return mask;
}

int cnt[LIM];
int magic(int n) /// O(n log max(a))
{
    memset(cnt, 0, sizeof(cnt[0]) * (n + 1));

    ll res = 0;
    for (int a = 1; a <= n; ++a) /// Check all cases of a
        res += cnt[getMask(a)]++;

    res %= MOD;
    return res;
}

int main() /// O(n log max(a))
{
    int n;
    cin >> n;
    sieve(n + 500);
    cout << magic(n);    
    return 0;
}

int solve(int n)
{
    memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));

    long long res = 0;
    for (int i = 1, j; i <= n; ++i) if (is_squarefree[i]) 
    {
        for (j = 1; i * j * j <= n; ++j)
            is_squarefree[i * j * j] = false;

        res += 1LL * (j - 1) * (j - 2) / 2;
    }

    res %= MOD;
    return res;
}

int solve(int n)
{
    memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));

    long long res = 0;
    for (int i = 1; i <= n; ++i) if (is_squarefree[i]) 
    {
        for (int j = 1; i * j * j <= n; ++j)
        {
            is_squarefree[i * j * j] = false;
            res += j - 1;
        }
    }

    res %= MOD;
    return res;
}

It just related to the sum of inversion of primes

Precalculating smallest prime factor for each $$$n$$$. And you can factorize any number in $$$O(log n)$$$

Each number will be calculated once

Though some numbers might be visited many times. It is actually armotized to constant

Let remind you about something called the sum of inversion of squares

Just ignore those squarefree. A normal loop through it is still linear.

$$$O\left(\underset{\text{squarefree } p \leq n}{\Large \Sigma} \left \lfloor {\Large \frac{n}{p^2}} \right \rfloor \right) = O\left(\underset{p \leq n}{\Large \Sigma} {\Large \frac{n}{p^2}} \right) = O\left(n \times {\Large \Sigma} {\Large \frac{1}{p^2}} \right) = O\left(n \times {\Large \frac{\pi^2}{6}} \right) = O(n)$$$

const int LIM = 1e7 + 17;
const int MOD = 1e9 + 7;

int fact[LIM + 1]; /// factorial:         fact[n] = n!
int invs[LIM + 1]; /// inverse modular:   invs[n] = n^(-1)
int tcaf[LIM + 1]; /// inverse factorial: tcaf[n] = (n!)^(-1)
void precal(int n = LIM)
{
    fact[0] = fact[1] = 1;
    invs[0] = invs[1] = 1;
    tcaf[0] = tcaf[1] = 1;
    for (int i = 2; i <= n; ++i)
    {
        fact[i] = (1LL * fact[i - 1] * i) % MOD;
        invs[i] = MOD - 1LL * (MOD / i) * invs[MOD % i] % MOD;
        tcaf[i] = (1LL * tcaf[i - 1] * invs[i]) % MOD;
    }
}

int nCk(int n, int k)
{
    k = min(k, n - k);
    if (k < 0) return 0;

    long long res = fact[n];
    res *= tcaf[k];         res %= MOD;
    res *= tcaf[n - k];     res %= MOD;
    return res;
}

int solve(int n)
{
    memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));

    long long res = 0;
    for (int i = 1, j; i <= n; ++i) if (is_squarefree[i]) 
    {
        for (j = 1; i * j * j <= n; ++j)
            is_squarefree[i * j * j] = false;

        res += nCk(j - 1, k);
    }

    res %= MOD;
    return res;
}

int solve(int n)
{
    memset(is_squarefree, true, sizeof(is_squarefree[0]) * (n + 1));
    long long res = 0;

    int t = sqrt(n);
    while (t * t < n) ++t;
    while (t * t > n) --t;
    for (int j = t; j > 1; --j)
    {
        for (int i = 1; i * j * j <= n; ++i)
        {
            if (!used[i * j * j])
            {
                used[i * j * j] = true;
                res += j - 1;
            }
        }
    }

    res %= MOD;
    return res;
}

#include <iostream>
#include <cstring>
#include <numeric>
#include <cmath>

using namespace std;

const int MOD = 1e9 + 7;
const int LIM = 1e7 + 17;
const int SQRT_LIM = ceil(sqrt(LIM) + 1) + 1;

int euler[SQRT_LIM];
void sieve_phi(int n)
{
    iota(euler, euler + n + 1, 0);
    for (int x = 2; x <= n; x++) if (euler[x] == x)
        for (int j = x; j <= n; j += x)
            euler[j] -= euler[j] / x;
}

int solve(int n)
{
    sieve_phi(ceil(sqrt(n) + 1) + 1);
    
    long long res = 0;
    for (int p = 2; p * p <= n; ++p)
        res += 1LL * euler[p] * (n / (p * p));

    res %= MOD;
    return res;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL);
    
    int n;
    cin >> n;
    cout << solve(n);
    return 0;
}

#include <iostream>
#include <cstring>
#include <numeric>
#include <vector>
#include <cmath>

using namespace std;

const int MOD = 1e9 + 7;

vector<int> lpf;
vector<int> prime;
vector<int> euler;
void linear_sieve_phi(int n)
{
    lpf.assign(n + 1, 0);
    euler.assign(n + 1, 1);
    for (int x = 2; x <= n; ++x)
    {
        if (lpf[x] == 0)
        {
            prime.push_back(lpf[x] = x);
            euler[x] = x - 1;                    
        }
        for (int i = 0; i < prime.size() && x * prime[i] <= n; ++i)
        {
            lpf[x * prime[i]] = prime[i];
            if (x % prime[i] == 0) {
                euler[x * prime[i]] = euler[x] * prime[i];    
                break;
            }
            euler[x * prime[i]] = euler[x] * euler[prime[i]];
        }
    }
}

int solve(int n)
{
    linear_sieve_phi(ceil(sqrt(n) + 1) + 1);
    
    long long res = 0;
    for (int p = 2; p * p <= n; ++p)
        res += 1LL * euler[p] * (n / (p * p));

    res %= MOD;
    return res;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL);
    
    int n;
    cin >> n;
    cout << solve(n);
    return 0;
}