大致是长剖+\(\rm dsu\ on\ tree\)的思想

先做一个转化，改为对于\(i\in[1,n-1]\)求出有多少个\(f(u,v)\)满足\(i|f(u,v)\)，这样我们最后再做一个反演就好了

既然我们要求有多少对\(f(u,v)\)是\(i\)或\(i\)的倍数，我们需要在长剖的时候快速合并两边的信息，这个信息长得非常别致，形如到当前节点距离为\(i\)或\(i\)的倍数的节点个数

轻儿子这边还好说，我们直接暴力调和级数处理一下即可，但是这样的信息从中儿子哪里却非常不好继承

考虑一下根号分治

如果这个\(i>\sqrt{n}\)，也就是直接暴力从重儿子那里暴力复杂度是\(\frac{n}{i}<\sqrt{n}\)，这样我们自然可以直接利用长剖时维护的数组得到

但是\(i\leq \sqrt{n}\)时，考虑维护一个数组可以快速得到这样的信息，但是这样的数组从重儿子那里没办法继承，看起来长剖好像行不太通

这个时候就需要用\(\rm dsu\)的思想了

考虑维护一个\(g[i][j]\)表示一个点子树内部点深度对\(i\)取模余数为\(j\)的点的个数，这样我们就没有必要考虑如何继承，像重儿子那样不清空直接拿来用就好了

利用这个数组我们很方便查询子树内部到当前节点距离为\(i\)或\(i\)的倍数的点的个数，只需要利用当前点的深度搞一搞就可以了

我们暴力轻儿子的时候也可以直接把更新\(g\)，更新一次复杂度是\(O(\sqrt{n})\)的

但一般来说，长剖都是先重儿子再轻儿子，但是\(\rm dsu\)却是先轻儿子再重儿子，这里为了保证\(g\)里维护的信息都来自于子树内部，我们必须得像\(\rm dsu\)一样先轻儿子之后处理重儿子

代码

#include
    
     #define re register#define LL long long#define min(a,b) ((a)<(b)?(a):(b))inline int read() {    char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();    while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;}const int maxn=2e5+5;struct E{int v,nxt;}e[maxn];int n,num,__,B;int head[maxn],son[maxn],len[maxn],dep[maxn],f[maxn],mu[maxn],p[maxn>>1];int h[maxn],top[maxn],dfn[maxn],g[320][320],b[maxn],c[maxn];LL ans[maxn],Ans[maxn];inline void add(int x,int y) {    e[++num].v=y;e[num].nxt=head[x];head[x]=num;}void dfs1(int x) {    for(re int i=head[x];i;i=e[i].nxt) {        dep[e[i].v]=dep[x]+1;dfs1(e[i].v);        if(len[e[i].v]>len[son[x]]) son[x]=e[i].v;    }    len[x]=len[son[x]]+1;}void dfs2(int x,int topf) {    top[x]=topf,dfn[x]=++__;    if(!son[x]) return;    dfs2(son[x],topf);    for(re int i=head[x];i;i=e[i].nxt)     if(son[x]!=e[i].v) dfs2(e[i].v,e[i].v);}inline void ins(int x,int v) {    for(re int i=1;i<=B;++i) g[i][x%i]+=v;}inline int calc(int x,int k) {    if(k<=B) return g[k][dep[x]%k];    int tot=0;    for(re int i=dfn[x]+k;i<=dfn[x]+len[x]-1;i+=k) tot+=h[i];    return tot;}void dfs(int x) {    for(re int i=head[x];i;i=e[i].nxt)     if(son[x]!=e[i].v) dfs(e[i].v);    if(son[x]) dfs(son[x]);    for(re int i=head[x];i;i=e[i].nxt) {        if(son[x]==e[i].v) continue;        int y=e[i].v;        for(re int j=1;j<=len[y];++j)            b[j]=h[dfn[y]+j-1];        for(re int j=1;j<=len[y];++j)            for(re int k=j;k<=len[y];k+=j) c[j]+=b[k];        for(re int j=1;j<=len[y];++j)            ans[j]+=1ll*c[j]*calc(x,j);        for(re int j=1;j<=len[y];++j) c[j]=b[j]=0;        for(re int j=1;j<=len[y];++j)            ins(j+dep[y]-1,h[dfn[y]+j-1]),h[dfn[x]+j]+=h[dfn[y]+j-1];    }    ins(dep[x],1);h[dfn[x]]++;    if(x==top[x]) {        for(re int i=1;i<=B;++i)            for(re int j=dep[x];j<=dep[x]+len[x]-1;++j)                g[i][j%i]=0;    }}int main() {    n=read();    for(re int x,i=2;i<=n;++i) x=read(),add(x,i);    dep[1]=1,dfs1(1),dfs2(1,1);f[1]=mu[1]=1;    B=std::ceil(std::sqrt(n));B=min(B,310);dfs(1);    for(re int i=2;i<=n;i++) {        if(!f[i]) p[++p[0]]=i,mu[i]=-1;        for(re int j=1;j<=p[0]&&p[j]*i<=n;++j) {            f[p[j]*i]=1;if(i%p[j]==0) break;            mu[p[j]*i]=-mu[i];        }    }    for(re int i=1;i<=n;i++)        for(re int j=i;j<=n;j+=i)            Ans[i]+=1ll*mu[j/i]*ans[j];    for(re int i=2;i<=n;++i)         c[1]++,c[dep[i]]--;    for(re int i=1;i<=n;i++) c[i]+=c[i-1];    for(re int i=1;i