NOIP训练 友好国度（点分治+容斥）

传送门
思路：
直接上点分治+容斥计算每个因数对应的贡献即可。
代码：

#include<bits/stdc++.h>
#define ri register int
using namespace std;
const int rlen=1<<18|1;
inline char gc(){
	static char buf[rlen],*ib,*ob;
	(ib==ob)&&(ob=(ib=buf)+fread(buf,1,rlen,stdin));
	return ib==ob?-1:*ib++;
}
inline int read(){
	int ans=0;
	char ch=gc();
	while(!isdigit(ch))ch=gc();
	while(isdigit(ch))ans=((ans<<2)+ans<<1)+(ch^48),ch=gc();
	return ans;
}
typedef long long ll;
const int N=1e5+5;
vector<int>e[N],fac[N],coe[N];
int n,a[N],lim=0,all,mx,siz[N],rt,Tim[N],cnt[N],mul[N],tim=0;
bool vis[N],chk[N];
ll ans=0;
inline int get(int x){
	if(mul[x])return mul[x];
	if(x<2)return mul[x]=x;
	int ret=1;
	for(ri i=0;i<fac[x].size();++i)ret*=fac[x][i];
	return mul[x]=ret;
}
inline void get_coef(int pos,int mul,int a){
	if(pos==fac[a].size()){if(~mul)coe[a].push_back(mul);return;}
	get_coef(pos+1,mul,a);
	get_coef(pos+1,-fac[a][pos]*mul,a);
}
inline void init(){
	for(ri i=2;i<=lim;++i)if(!vis[i])for(ri j=i;j<=lim;j+=i)vis[j]=1,fac[j].push_back(i);
	for(ri i=1,x;i<=n;++i)vis[i]=0,a[i]=get(a[i]);
}
void getroot(int p,int fa){
	siz[p]=1;
	int ms=0;
	for(ri i=0,v;i<e[p].size();++i){
		if((v=e[p][i])==fa||vis[v])continue;
		getroot(v,p),siz[p]+=siz[v],ms=max(ms,siz[v]);
	}
	ms=max(ms,all-siz[p]);
	if(ms<=mx)rt=p,mx=ms;
}
inline int query(int x){return Tim[x]==tim?cnt[x]:0;}
inline int ask(int x){
	int ret=0;
	if(!chk[x])get_coef(0,-1,x),chk[x]=1;
	for(ri i=0;i<coe[x].size();++i)ret+=coe[x][i]/abs(coe[x][i])*query(abs(coe[x][i]));
	return ret;
}
inline void update(int x){Tim[x]==tim?++cnt[x]:cnt[x]=1,Tim[x]=tim;}
inline void change(int x){
	if(!chk[x])get_coef(0,-1,x),chk[x]=1;
	for(ri i=0;i<coe[x].size();++i)update(abs(coe[x][i]));
}
inline int gcd(int a,int b){int t;while(b){t=a,a=b,b=t-t/a*a;}return a;}
void dfs(int p,int fa,int g,int coef){
	g=gcd(g,a[p]),ans+=coef*ask(g),change(g),siz[p]=1;
	for(ri i=0,v;i<e[p].size();++i){
		if(vis[v=e[p][i]]||v==fa)continue;
		dfs(v,p,g,coef),siz[p]+=siz[v];
	}
}
inline void solve(int p){
	vis[p]=1;
	++tim;
	change(a[p]);
	for(ri i=0,v;i<e[p].size();++i)if(!vis[v=e[p][i]])dfs(v,p,a[p],1);
	for(ri i=0,v;i<e[p].size();++i){
		if(vis[v=e[p][i]])continue;
		++tim,dfs(v,p,a[p],-1),mx=all=siz[v],getroot(v,p),solve(rt);
	}
}
int main(){
	n=read();
	for(ri i=1,u,v;i<n;++i)u=read(),v=read(),e[u].push_back(v),e[v].push_back(u);
	for(ri i=1;i<=n;++i)a[i]=read(),lim=max(lim,a[i]);
	init();
	rt=0,mx=all=n,getroot(1,0),solve(rt);
	cout<<(ll)n*(n-1)/2-ans;
	return 0;
}