Consider the following function \(\tilde{\mu}(n)\) that takes a positive integer \(n\) and returns either of \(\{ -1,0,1\}\). (This function is based on Möbius function. For details, please refer to the appendix.)

\[ \tilde{\mu}(n)= \begin{cases} 1 & (\text{if }n\text{ is expressed by odd number of distinct primes }(\text{like }n=2,3,30) ) \\ -1 & (\text{if }n\text{ is expressed by even (and positive) number of distinct primes }(\text{like }n=6,10,210) ) \\ 0 & (\text{otherwise }(\text{like }n=1,4,8,9,12) ) \end{cases} \]

This function satisfies the following property.

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

Here, \(\displaystyle\sum_{d|n} \mu(d)\) denotes the sum of \(\mu(d)\) over all divisors \(d\) of \(n\).
Next, for the integers \(a\) and \(b\) such that \(1\leq a\leq N\) and \(1\leq b\leq N\), and for the pairs of integers such that \(1\leq i\leq j\leq N\), consider the following function \(f(a,b;i,j)\).

\[ f(a,b;i,j)= \begin{cases} 1 & ( \text{if }(i\text{ is a multiple of }a)\text{ and }(p_i\text{ is a multiple of }b)\text{ and }(j\text{ is a multiple of }a)\text{ and }(p_j\text{ is a multiple of }b) ) \\ 0 & ( \text{otherwise} ) \end{cases} \]

Then, consider

\[ 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). \]

For each \((i, j)\), let \(GCD(i,j)=g\) and \(GCD(p_i,p_j)=g'\), then

\[ \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), \]

and the right hand side is \(1\) if both \(g\) and \(g'\) are greater than or equal to \(2\), or otherwise \(0\), so \(S\), the sum of this expression over all \(1\leq i\leq j\leq N\), is equal to the desired value.

On the other hand, let \(num(a,b)\) be the number of integers \(i\) such that \(1\leq i\leq N\) and (\(i\) is a multiple of \(a\)) and (\(p_i\) is a multiple of \(b\)), then \(\displaystyle\sum_{1\leq i\leq j\leq N} f(a,b;i,j)=\frac{num(a,b)(num(a,b)+1)}{2}\), so it holds that

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

If you try to compute this naively, it will cost \(O(N^2)\) time, but actually, for each \(a\), it is sufficient to compute the value for only \(b\) such that for some multiple \(i\) of \(a\), \(b\) is a divisor of \(p_i\). Moreover, it is sufficient to check only those such that \(\tilde{\mu}(b)\neq 0\), and since \(2\cdot 3\cdot 5\cdot 7\cdot 11\cdot 13\cdot 17=510510>200000\), for each \(p_i\) the number of such \(b\) is at most \(2^6-1=63\). Therefore, the number of computations is at most about \(\displaystyle\sum_{a=1}^N \left( 63\cdot \frac{N}{a}\right)\simeq 63N\log N\simeq 2\times 10^8\), and since the constant factor is not really heavy, so it will fit in the time limit. Therefore, the problem has been solved.

Appendix

The Möbius function \(\mu(n)\) is defined as

\[ \tilde{\mu}(n)= \begin{cases} 1 & (\text{if }n\text{ is expressed by even (possibly zero) number of distinct primes }(\text{like }n=1,6,10,210) ) \\ -1 & (\text{if }n\text{ is expressed by odd number of distinct primes }(\text{like }n=2,3,30) ) \\ 0 & (\text{otherwise }(\text{like }n=4,8,9,12) ) \end{cases} \]

and it holds that

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

Möbius function often appears in the problem that treats inclusion-exclusion of divisors. This time, we want to inspect those elements with the greatest common divisor greater than \(1\), so it is convenient that the function takes the value \(1\) for \(n\geq 2\). That is why we defined

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

so that \(\displaystyle\sum_{d|n} \tilde{\mu}(d)=1-\sum_{d|n}\mu(d)\).

Also, since \(\displaystyle\sum_{a=1}^N \sum_{b=1}^N g(a,b)^2\sim O(N^2\log N)\), as long as the computation is done in 64-bit integers, the overflow will never happen.

Sample code in 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;
}