前置知识：莫比乌斯反演。

首先，考虑对于一组 $(i,j)$，令 $g=\gcd(i,j),g'=\gcd(p_i,p_j)$。

直接用莫比乌斯反演一下。

令 $f(a,b,i,j)$ 表示 $(i,j)$ 是否满足 $g=a,g'=b$，可以是 $1$，不可以是 $0$。

然后发现，对于一组 $(i,j)$，满足条件的 $a,b$ 实际上就是我们在上面那个式子枚举的 $a,b$。

这是一个比较复杂的式子，我们考虑化简。

设 $s(a,b)$ 表示对于一组 $a,b$，有多少个 $i$ 满足 $i$ 是 $a$ 的倍数，$p_i$ 是 $b$ 的倍数。

可以得到

然后去化简原式，得到

然后就可以发现，只有满足 $i$ 是 $a$ 的倍数的位置才会对答案产生贡献，这样的位置只有 $\frac{n}{a}$ 个，此时 $b$ 要满足是这些位置的 $p_i$ 的因数。

由于 $p_i\le 2\times 10^5$，对于一个位置，它的因数最多只有 $63$ 个，先预处理这些位置，直接枚举就好了。

时间复杂度 $O(63n\log n)$。

#include<bits/stdc++.h>
using namespace std;
const int N=2e5+5;
int n,p[N],mu[N],b[N],P[N],cnt,s[N],vis[N];
vector<int> cana,vec[N],canb;
long long ans;
void work(int maxn){//预处理莫比乌斯函数
	mu[1]=0;
	for(int i=2;i<=maxn;i++){
		if(!b[i])P[++cnt]=i,mu[i]=1;
		for(int j=1;j<=cnt&&P[j]*i<=maxn;j++){
			b[P[j]*i]=1;
			if(i%P[j]==0){
				mu[i*P[j]]=0;
				break;
			}
			mu[i*P[j]]=-mu[i];
		}
		if(mu[i]){//如果值为零，对答案显然没有贡献，所以只考虑不为零的
			cana.push_back(i);
			for(int j=i;j<=maxn;j+=i)vec[j].push_back(i);
		}
	}
}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;i++)scanf("%d",&p[i]);
	work(n);
	for(int a:cana){
		canb.clear();
		for(int i=a;i<=n;i+=a){
			for(int b:vec[p[i]]){
				s[b]++;
				if(!vis[b])vis[b]=1,canb.push_back(b);
			}
		}
		for(int b:canb){
			ans=ans+1ll*mu[a]*mu[b]*s[b]*(s[b]+1)/2;
			s[b]=0;
			vis[b]=0;
		}
	}
	printf("%lld",ans);
} 

题意

给定一个长度为 $n$ 的排列 $p$，求有多少 $1\le i\le j\le n$ 满足 $\gcd(i,j)=1,\gcd(p_i,p_j)=1$。

分析

看到 $\gcd$ 计数题，我们考虑莫比乌斯反演。 然后就开始推：

其中 $f(a,b)=\sum_{1\le i\le n}[a|i][b|p_i]$。

然后，由于在本题范围内一个数的 $\mu$ 值不为 $0$ 的约数有最多有 $63$ 个，直接枚举 $a$，然后枚举 $i$ 和 $b$ 并统计每个 $b$ 的 $num$，最后枚举 $b$ 统计答案。

#include<bits/stdc++.h>
using namespace std;
#define ll long long
int n;
int per[200005];
vector<int> d[200005],vec,v;
int mu[200005];
int check[200005];
int num[200005];
int prime[200005],tot;
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;i++)scanf("%d",&per[i]);
	for(int i=2;i<=n;i++){
		if(!check[i]){
			mu[i]=-1;
			prime[++tot]=i;
		}
		for(int j=1;j<=tot&&prime[j]<=n/i;j++){
			check[i*prime[j]]=1;
			if(i%prime[j]==0){
				mu[i*prime[j]]=0;
				break;
			}else mu[i*prime[j]]=-mu[i];
		}
		if(mu[i]){
			vec.push_back(i);
			for(int j=i;j<=n;j+=i)d[j].push_back(i);
		}
	}
	ll ans=0;
	for(int a:vec){
		v.clear();
		for(int i=a;i<=n;i+=a)
			for(int b:d[per[i]]){
				if(!num[b])v.push_back(b);
				num[b]++;
			}
		for(int b:v){
			ans+=1ll*mu[a]*mu[b]*num[b]*(num[b]+1)/2;
//			printf("%d %d %lld\n",a,b,ans);
			num[b]=0;
		}
	}
	printf("%lld\n",ans);
	return 0;
}