Link
考虑计算出每个点除最后一步以外的期望到达次数，很显然初始同色点的这个值是相等的。
设(f_{i,0/1})表示在有(i)个(1)时，颜色为(0/1)的点除最后一步以外的期望到达次数。
那么可以列出转移方程：
[egin{aligned}
f_{i,0}&=frac inf_{i-1,0}+frac{n-i-1}nf_{i+1,0}+frac1nf_{i+1,1}+frac1n\
f_{i,1}&=frac{n-i}nf_{i+1,1}+frac{i-1}nf_{i-1,1}+frac1nf_{i-1,0}+frac1n
end{aligned}
]
利用第二个方程计算出(f_{i+1,1})，再利用第一个方程计算出(f_{i+1,0})即可完成递推。
那么设(f_{1,0}=x,f_{1,1}=y)即可将所有的(f_{i,o})都用(ax+by+c)表示。
最后联立两个方程解出(x,y)即可计算得到答案。
#include<cctype>
#include<cstdio>
#include<vector>
using i64=long long;
const int N=100007,P=1000000007;
char col[N];int n,fa[N],size[N];i64 inv[N],ans[N];
std::vector<int>e[N];
struct node
{
i64 a,b,c;
node(int A=0,int B=0,int C=0):a(A),b(B),c(C){}
void mod(){a+=a>>63&P,b+=b>>63&P,c+=c>>63&P;}
}f[N][2];
node operator+(const node a,const node b){return node(a.a+b.a,a.b+b.b,a.c+b.c);}
node operator-(const node a,const node b){return node(a.a-b.a,a.b-b.b,a.c-b.c);}
node operator*(const node a,const i64 b){return node(a.a*b%P,a.b*b%P,a.c*b%P);}
node operator+(const node a,const i64 b){return node(a.a,a.b,a.c+b);}
node operator-(const node a,const i64 b){return node(a.a,a.b,a.c-b);}
int read(){int x=0,c=getchar();while(isspace(c))c=getchar();while(isdigit(c))(x*=10)+=c&15,c=getchar();return x;}
i64 pow(i64 a,i64 b){a%=P;i64 r=1;for(;b;b>>=1,a=a*a%P)if(b&1)r=r*a%P;return r;}
void dfs1(int u)
{
size[u]=1;
for(int v:e[u]) if(v^fa[u]) dfs1(v),size[u]+=size[v],ans[u]+=ans[v]+size[v];
}
void dfs2(int u)
{
if(fa[u]) ans[u]=ans[fa[u]]+n-2*size[u];
for(int v:e[u]) dfs2(v);
}
int main()
{
int cnt=0;
n=read(),scanf("%s",col+1),inv[0]=inv[1]=f[1][0].a=f[1][1].b=1;
for(int i=1;i<=n;++i) cnt+=col[i]=='1';
for(int i=2;i<=n;++i) e[fa[i]=read()].push_back(i),inv[i]=(P-P/i)*inv[P%i]%P;
dfs1(1),dfs2(1);
for(int i=1;i<n;++i) f[i+1][1]=(f[i][1]*n-f[i-1][1]*(i-1)-f[i-1][0]-(i!=1))*inv[n-i],f[i+1][0]=(f[i][0]*n-f[i-1][0]*i-f[i+1][1]-1)*inv[n-i-1],f[i+1][0].mod(),f[i+1][1].mod();
node a=f[n-1][1]*n-f[n-2][1]*(n-2)-f[n-2][0]-1,b=f[n-1][0]*n-f[n-2][0]*(n-1);
i64 x=(b.b*a.c-b.c*a.b)%P*pow(a.b*b.a-a.a*b.b,P-2)%P,y=(b.a*a.c-b.c*a.a)%P*pow(a.a*b.b-a.b*b.a,P-2)%P,coef[2],Ans=0;
for(int i=0;i<2;++i) coef[i]=(f[cnt][i].a*x+f[cnt][i].b*y+f[cnt][i].c+inv[n])%P;
for(int i=1;i<=n;++i) (Ans+=ans[i]%P*coef[col[i]&15])%=P;
printf("%lld",(Ans+P)*inv[n]%P);
}