この問題を \(\Theta(\sqrt{N})\) 時間で解く方法について解説します。

答えは \(\displaystyle2\sum_{d=2}^{\lfloor\sqrt{N}\rfloor}\left\lfloor\dfrac N{d^2}\right\rfloor\varphi(d)+N\) です。ただし、\(\varphi(n)\) は \(n\) と互いに素な \(n\) 以下の正の整数の個数です。
\(\varphi(2),\varphi(3),\ldots,\varphi(\lfloor\sqrt{N}\rfloor)\) は線形篩を用いて \(\Theta(\sqrt{N})\) 時間で求められるため、全体でも \(\Theta(\sqrt{N})\) 時間でこの問題が解けました。以下でこの式の導出と具体的な求め方を説明します。

導出

数える対象である正整数の組 \((i,j)\) に対して、\(g\coloneqq\gcd(i,j),i^\prime\coloneqq\dfrac{i}{g},j^\prime\coloneqq\dfrac{j}{g}\) と定めます。 このとき \(i^\prime,j^\prime\) はいずれも平方数で、互いに素です。
逆に、\(1\leq ga^2\leq N, 1\leq gb^2\leq N, \gcd(a,b)=1\) を満たす正整数の \(3\) つ組 \((g,a,b)\) を任意に取ったとき \(i=ga^2,j=gb^2\) は所望の条件を満たします。

よって、この問題は \(1\leq ga^2\leq N, 1\leq gb^2\leq N, \gcd(a,b)=1\) を満たす正整数の \(3\) つ組 \((g,a,b)\) の個数を求めることに帰着されました。

\(a=b=1\) である \(N\) 個は別に考え、\(a\lt b\) の場合の数を倍することを考えます。

\(b\) を固定すると、

• \(a\) は \(b\) と互いに素な \(b\) 以下の正の整数であるため、これは \(\varphi(b)\) 個あります。
• \(g\) は \(\dfrac{N}{b^2}\) 以下の正の整数を自由に動けるため、これは \(\left\lfloor\dfrac{N}{b^2}\right\rfloor\) 個あります。

\(a,g\) は（それぞれが条件を満たせば）自由に選べるため、\(b\) を固定したとき条件を満たす \(3\) つ組の個数は \(\left\lfloor\dfrac{N}{b^2}\right\rfloor\varphi(b)\) 個です。

（\(a\lt b\) も含めた）条件を満たす \(3\) つ組において、 \(b\) は \(2\leq b\leq\sqrt{N}\) を満たします。 なので、（\(a\lt b\) も含めた）条件を満たす \(3\) つ組の個数は \(\displaystyle\sum_{d=2}^{\lfloor\sqrt{N}\rfloor}\left\lfloor\dfrac N{d^2}\right\rfloor\varphi(d)\) 個です。

対称性より、\(a\gt b\) を満たす \(3\) つ組は \(a\lt b\) を満たす \(3\) つ組と同じだけあります。
\(a=b\) を満たす \(3\) つ組は \(N\) 個あったので、\(\displaystyle2\sum_{d=2}^{\lfloor\sqrt{N}\rfloor}\left\lfloor\dfrac N{d^2}\right\rfloor\varphi(d)+N\) が答えであることがわかりました。

こちらの解説における \(\displaystyle\sum_d|S_d|d^2=\sum_d|T_d|(2d-1)\) を求める問題だと考えたのち、除原理の操作の逆を \(2d-1\) に対して行って変形していくことでも同じ式を導出することができます。

求め方

線形篩の詳細はこちらの解説にあるリンクを参考にしてください。

線形篩は、\(\dfrac{n}{\operatorname{lpf}_n}\rightarrow n\) なる遷移のみを考えることで \(i=2,3,\ldots,X\) についての \(\operatorname{lpf}_i\coloneqq i\) の最小素因数 を時間計算量 \(\Theta(X)\) で全て求めるアルゴリズムです。

いま、線形篩を用いて \(\operatorname{lpf}_2,\operatorname{lpf}_3,\ldots,\operatorname{lpf}_{\lfloor\sqrt{N}\rfloor}\) が求められているとします。

このとき、次の漸化式にしたがって計算した \(\operatorname{dp}[i]\) が \(\varphi(i)\) を与えます。

\[\begin{aligned} \operatorname{dp}[1]&=1 \\ \operatorname{dp}[i]&=\left\{\begin{array}{ll}i-1&(\operatorname{lpf}_i=i)\\\operatorname{dp}\left[\dfrac{i}{\operatorname{lpf}_i}\right]\times\operatorname{lpf}_i&(\operatorname{lpf}_i=\operatorname{lpf}_{i/\operatorname{lpf}_i})\\\\\operatorname{dp}\left[\dfrac{i}{\operatorname{lpf}_i}\right]\times(\operatorname{lpf}_i-1)&(\text{otherwise.})\end{array}\right. \end{aligned}\]

あとは、この \(\operatorname{dp}[i]\) を用いて \(\displaystyle2\sum_{d=2}^{\lfloor\sqrt{N}\rfloor}\left\lfloor\dfrac N{d^2}\right\rfloor\varphi(d)+N\) を計算すればよいです。

実装は以下のようになります。(実際の提出)

#include <bits/stdc++.h>
using namespace std;

int main(){
    int N;
    cin >> N;
    int sqrt_N = sqrt(N);

    // 線形篩で前計算 Θ(√N)時間
    vector<int> primes, least_prime_factor(sqrt_N + 1);
    for(int i = 2; i <= sqrt_N; ++i){
        if(!least_prime_factor[i]){
            primes.emplace_back(i);
            least_prime_factor[i] = i;
        }
        for(const auto p : primes){
            if(p * i > sqrt_N || p > least_prime_factor[i])break;
            least_prime_factor[p * i] = p;
        }
    }

    // φ(2), φ(3), ..., φ(Floor(√N)) を線形篩を用いて計算 Θ(√N)時間
    vector<int> totient(sqrt_N + 1);
    for(int i = 2; i <= sqrt_N; ++i){
        const auto p = least_prime_factor[i];
        totient[i] = i == p ? i - 1
                   : least_prime_factor[i / p] == p ? p * totient[i / p]
                   : (p - 1) * totient[i / p];
    }

    // 答えは 2 * ∑ Floor(N/i²) φ(i) + N
    int ans = 0;
    for(int i = 2; i <= sqrt_N; ++i)ans += N / i / i * totient[i];
    cout << 2 * ans + N << endl;

    return 0;
}