Link
考虑计算出每个点除最后一步以外的期望到达次数，很显然初始同色点的这个值是相等的。
设\(f_{i,0/1}\)表示在有\(i\)个\(1\)时，颜色为\(0/1\)的点除最后一步以外的期望到达次数。
那么可以列出转移方程：

\[\begin{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);
}