正整数 \(n\) に対して \(\{ -1,0,1\}\) のいずれかを返す次のような関数 \(\tilde{\mu}(n)\) を考えます。（これはメビウス関数を元に作られています。詳しくは下の補足を見てください。）

\[ \tilde{\mu}(n)= \begin{cases} 1 & (nが相異なる奇数個の素数の積で表されるとき (n=2,3,30 など) ) \\ -1 & (nが相異なる偶数個 (0 を含まない) の素数の積で表されるとき (n=6,10,210 など) ) \\ 0 & (それ以外のとき (n=1,4,8,9,12 など) ) \end{cases} \]

そしてこれは次の性質をみたしています。

\[ \sum_{d|n} \tilde{\mu}(d)= \begin{cases} 0 & (n=1 のとき) \\ 1 & (n\geq 2 のとき) \end{cases} \]

ただし、\(\displaystyle\sum_{d|n} \mu(d)\) は \(n\) の約数である \(d\) について \(\mu(d)\) を足し合わせることを表しています。
次に、\(1\leq a\leq N\) , \(1\leq b\leq N\) をみたす整数 \(a,b\) および \(1\leq i\leq j\leq N\) をみたす整数の組 \((i,j)\) に対して、次の関数 \(f(a,b;i,j)\) を考えます。

\[ f(a,b;i,j)= \begin{cases} 1 & ( (i が a の倍数) かつ (p_i が b の倍数) かつ (j が a の倍数) かつ (p_j が b の倍数) のとき ) \\ 0 & ( それ以外のとき) \end{cases} \]

このとき、

\[ S=\sum_{1\leq i\leq j\leq N} \sum_{a=1}^N \sum_{b=1}^N \tilde{\mu}(a) \tilde{\mu}(b)f(a,b;i,j) \]

を考えます。各 \((i,j)\) に対して、\(GCD(i,j)=g\), \(GCD(p_i,p_j)=g'\) とすると、

\[ \sum_{a=1}^N \sum_{b=1}^N \tilde{\mu}(a) \tilde{\mu}(b)f(a,b;i,j) =\sum_{a|g}\sum_{b|g'} \tilde{\mu}(a) \tilde{\mu}(b) \]

であり、この右辺は \(g\) と \(g'\) がともに \(2\) 以上のとき \(1\) 、そうでないとき \(0\) となるためこれを \(1\leq i\leq j\leq N\) で足し合わせた \(S\) は求めるべき値と一致します。

一方で、 \(num(a,b)\) を \(1\leq i\leq N\) かつ ( \(i\)が \(a\) の倍数 ) かつ ( \(p_i\)が \(b\) の倍数 ) をみたす整数 \(i\) の個数とすると、 \(\displaystyle\sum_{1\leq i\leq j\leq N} f(a,b;i,j)=\frac{num(a,b)(num(a,b)+1)}{2}\) より、

\[ S=\sum_{a=1}^N \sum_{b=1}^N \tilde{\mu}(a) \tilde{\mu}(b) \frac{num(a,b)(num(a,b)+1)}{2} \]

と書けます。

これを愚直に計算しようとすると \(O(N^2)\) かかってしまいますが、実際には各 \(a\) について、ある\(a\)の倍数 \(i\) について \(p_i\) の約数であるような \(b\) についてだけ \(num(a,b)\) を計算すればよいです。 さらに \(\tilde{\mu}(b)\neq 0\) であるものについてだけ調べればよく、\(2\cdot 3\cdot 5\cdot 7\cdot 11\cdot 13\cdot 17=510510>200000\) より、各 \(p_i\)について そのような \(b\) の個数は高々 \(2^6-1=63\) 個です。よって計算回数は高々 \(\displaystyle\sum_{a=1}^N \left( 63\cdot \frac{N}{a}\right)\simeq 63N\log N\simeq 2\times 10^8\) 回程度であり、定数倍も重くないので十分間に合います。 よって、この問題を解く事が出来ました。

補足

メビウス関数 \(\mu(n)\) は

\[ \mu(n)= \begin{cases} 1 & (nが相異なる偶数個 (0 を含む) の素数の積で表されるとき (n=1,6,10,210 など) ) \\ -1 & (nが相異なる奇数個の素数の積で表されるとき (n=2,3,30 など) ) \\ 0 & (それ以外のとき (n=4,8,9,12 など) ) \end{cases} \]

で定義され、

\[ \sum_{d|n}\mu(d)= \begin{cases} 1 & (n=1 のとき) \\ 0 & (n\geq 2 のとき) \end{cases} \]

が成り立ちます。

メビウス関数は約数包除を使う問題でよく出てますが、 今回は、ともに最大公約数が \(1\) でないものを調べたいため \(n\geq 2\) で \(1\) となる方が都合がよく、

\[ \tilde{\mu}(n)= \begin{cases} 1-\mu (n) & (n=1のとき) \\ -\mu(n) & (n\geq 2のとき) \end{cases} \]

と定めることで、\(\displaystyle\sum_{d|n} \tilde{\mu}(d)=1-\sum_{d|n}\mu(d)\) となります。

また、\(\displaystyle\sum_{a=1}^N \sum_{b=1}^N g(a,b)^2\sim O(N^2\log N)\) であるため、 \(64\) bit 整数型で計算を行っている限り、計算中にオーバーフローをすることもありません。

c++による実装例:

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

#define N 200010

int main() {
	int n;
	int p[N];

	bool pr[N];
	int mu[N];
	vector<int>m;
	vector<int>d[N];

	bool used[N];
	vector<int>cand;
	long long num[N];

	long long ans = 0;

	cin >> n;
	for (int i = 0; i < n; i++)cin >> p[i + 1];

	for (int i = 1; i <= n; i++) {
		pr[i] = true;
		mu[i] = -1;
	}
	pr[1] = false;
	mu[1] = 0;

	for (int i = 2; i <= n; i++) {
		if (pr[i]) {
			mu[i] = -mu[i];
			for (int j = 2; (i*j) <= n; j++) {
				pr[i*j] = false;
				if (j%i == 0)mu[i*j] = 0;
				else mu[i*j] = -mu[i*j];
			}
		}
		if (mu[i] != 0) {
			m.push_back(i);
			for (int j = i; j <= n; j += i)d[j].push_back(i);
		}
	}

	for (int i = 1; i <= n; i++) {
		used[i] = false;
		num[i] = 0;
	}
	for (auto &a : m) {
		for (int i = a; i <= n; i += a) {
			for (auto &j : d[p[i]]) {
				num[j]++;
				if (!used[j]) {
					used[j] = true;
					cand.push_back(j);
				}
			}
		}
		for (auto &b : cand) {
			ans += mu[a] * mu[b] * ((num[b] * (num[b] + 1)) / 2);
			num[b] = 0;
			used[b] = false;
		}
		cand.clear();
	}

	cout << ans << endl;


	return 0;
}