Author: Kamonasish Roy
Tester : Emon Thakur
Editorial :
Key observations
Let d = |i − j|.
The equation gcd(i, j) = d implies that d divides both i and j, therefore

 i = k·d, j = (k + 1)·d for some integer k ≥ 1.

Thus, for a fixed difference d a valid k must satisfy

 (k + 1)·d ≤ N ⇒ k ≤ ⌊N/d⌋ − 1.

Consequently the number of harmonious pairs for this d is

 ⌊N/d⌋ − 1.

Summing over all d from 1 to N,

 Answer(N) = Σ₍d = 1…N₎ (⌊N/d⌋ − 1)
= Σ₍d=1…N₎ ⌊N/d⌋ − N.

Define
H(N) = Σ₍d=1…N₎ ⌊N/d⌋.
Then
Answer(N) = H(N) − N.

Computing H(N) in O(√N)
Naively evaluating H(N) costs O(N), which is impossible for N = \(10^9\).
We exploit the classical “floor‑sum” trick.

Let m = ⌊√N⌋.

Direct summation for small divisors
For k = 1..m we can compute ⌊N/k⌋ directly.
Call partial = Σ₍k=1…m₎ ⌊N/k⌋.

Symmetry for large divisors
Every integer value v = ⌊N/k⌋ with v ≤ m appears the same
number of times on the other side of the square root.
A proven identity (derivable from grouping equal quotients) is

     H(N) = 2·partial − m².

Therefore

 Answer(N) = (2·partial − m²) − N.

Because partial needs only O(m) iterations, the total complexity per test
case is O(√N) time and O(1) memory.

Code (C++) :

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

ll harmonious_pairs(ll N) {
    ll m = (ll)std::sqrt((long double)N);
    ll partial = 0;
    for (ll k = 1; k <= m; ++k) partial += N / k;
    ll H = (partial << 1) - m * m;   // 2*partial - m^2
    return H - N;                    // subtract N per formula
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int T; cin >> T;
    while (T--) {
        ll N; cin >> N;
        cout << harmonious_pairs(N) << '\n';
    }
    return 0;
}