Denote $cnt_{x}$ = the number of occurrences of $x$, $h$ = the maximum of $a_i$, there we get
$$ \sum_{1\leqslant i\leqslant n}\sum_{1\leqslant j\leqslant n}[a_i,a_j] \\ \begin{aligned} &=\sum_{1\leqslant d\leqslant h}\sum_{1\leqslant i\leqslant h}\sum_{1\leqslant j\leqslant h}[\left(i,j\right)=d]\times\frac{i\times j\times cnt_i\times cnt_j}{d} \\ &=\sum_{1\leqslant d\leqslant h}d\sum_{1\leqslant i\leqslant\lfloor\frac{h}{d}\rfloor}\sum_{1\leqslant j\leqslant\lfloor\frac{h}{d}\rfloor}\sum_{k\mid\left(i,j\right)}\mu\left(k\right)\times i\times j\times cnt_{id}\times cnt_{jd} \\ &=\sum_{1\leqslant d\leqslant h}d\sum_{1\leqslant k\leqslant\lfloor\frac{h}{d}\rfloor}\mu\left(k\right)\times k^2\sum_{1\leqslant i\leqslant\lfloor\frac{h}{dk}\rfloor}\sum_{1\leqslant j\leqslant\lfloor\frac{h}{dk}\rfloor}i\times j\times cnt_{idk}\times cnt_{jdk} \\ &=\sum_{1\leqslant T\leqslant h}T\sum_{k\mid T}\mu\left(k\right)\times k\sum_{1\leqslant i\leqslant \lfloor\frac{h}{T}\rfloor}i\times cnt_{iT}\sum_{1\leqslant j\leqslant \lfloor\frac{h}{T}\rfloor}j\times cnt_{jT} \\ \end{aligned} $$
Denote $\displaystyle f(T)=\sum_{1\leqslant i\leqslant h/T}i\times cnt_{iT}$, $\displaystyle g(T)=\sum_{1\leqslant i\leqslant h/T}i\times cnt_{iT}\times f(T)$
$$ \sum_{1\leqslant T\leqslant h}T\sum_{k\mid T}\mu\left(k\right)\times k\sum_{1\leqslant i\leqslant \lfloor\frac{h}{T}\rfloor}i\times cnt_{iT}\sum_{1\leqslant j\leqslant \lfloor\frac{h}{T}\rfloor}j\times cnt_{jT} \\ \begin{aligned} &=\sum_{1\leqslant T\leqslant h}T\sum_{k\mid T}\mu\left(k\right)\times k\times g\left(T\right) \end{aligned} $$
Denote $\displaystyle z(T)=T\sum_{k\mid T}\mu(k)\times k$, the answer would be $\displaystyle\sum_{1\leqslant i\leqslant h}z(i)\times g(i)$.
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define cmin(x, y) x = min(x, y)
#define cmax(x, y) x = max(x, y)
#define fors(i, l, r, ...) for(int i = (l), REP##i = (r), ##__VA_ARGS__; i <= REP##i; ++i)
#define dfors(i, r, l, ...) for(int i = (r), REP##i = (l), ##__VA_ARGS__; i >= REP##i; --i)
signed main() {
ios::sync_with_stdio(0), cin.tie(0);
int n;
cin >> n;
vector<int> a(n);
for(int& x : a) cin >> x;
const int h = *max_element(a.begin(), a.end());
vector<int> cnt(h+1);
for(const int x : a) cnt[x]++;
const auto z = [](const int n) {
vector<int> mu(n+1), prime, not_prime(n+1);
vector<ll> z(n+1);
mu[1] = 1;
fors(i, 2, n) {
if(not_prime[i] == 0) prime.emplace_back(i),mu[i] = -1;
for(const int p : prime) {
if(i > n/p) break;
not_prime[i*p] = 1;
if(i%p == 0) break;
mu[i*p] = -mu[i];
}
}
fors(d, 1, n) {
for(int T = d; T <= n; T += d) z[T] += ll(mu[d])*d;
}
return z;
}(h);
ll ans = 0;
fors(i, 1, h) {
ll tmp = 0;
fors(j, 1, h/i) tmp += ll(cnt[i*j])*j;
ans += tmp*tmp*z[i]*i;
}
cout << ans << "\n";
return 0;
}