#if defined(_USE_PCH_)
#include "pch.hpp"
#else
#include <bits/stdc++.h>
#endif
#define RNG(V_, A_, B_, ...) for(int V_=(A_), V_##_END=(B_) __VA_OPT__(,) __VA_ARGS__; V_<=V_##_END; V_++)
#define IRNG(V_, A_, B_, ...) for(int V_=(A_), V_##_END=(B_) __VA_OPT__(,) __VA_ARGS__; V_>=V_##_END; V_--)
#ifdef _WIN32
#define long int64_t
#endif
#define fir first
#define sec second
#define _UN using namespace
using namespace std;
typedef pair<int,int> pii;
typedef __int128 i128;
const int MAXN=5e5+10; const long MOD=1e9+7;
int n,k,fa[MAXN],A[MAXN]; vector<int> G[MAXN];
char tmp[MAXN];
void addto(long& x,long y){ x=(x+y)%MOD; }
long qpow(long x,int p){
    long r=1;
    for(; p; p>>=1,x=x*x%MOD){
        if(p&1) r=r*x%MOD;
    }
    return r;
}
struct X{
    long a,b,c; // ax+by+c
    X operator+(X o)const{ return {(a+o.a)%MOD,(b+o.b)%MOD,(c+o.c)%MOD}; }
    X operator-(X o)const{ return {(a-o.a+MOD)%MOD,(b-o.b+MOD)%MOD,(c-o.c+MOD)%MOD}; }
    X operator*(long x)const{ return {(a*x)%MOD,(b*x)%MOD,(c*x)%MOD}; }
} F[MAXN][2]; long H[MAXN][2];
long D[MAXN],D2[MAXN]; int siz[MAXN];
void dfs1(int u){
    siz[u]=1;
    for(int v:G[u]) dfs1(v),D2[u]+=D2[v]+siz[v];
}
void dfs2(int u,long d){
    D[u]=(D2[u]+d)*qpow(n,MOD-2);
    for(int v:G[u]) dfs2(v,d+(D2[u]-(D2[v]+siz[v]))+(n-siz[v]));
}
int main(){
#if defined(_LOCAL_)
    freopen("in","r",stdin);
//  freopen("out","w",stdout);
//  freopen("/dev/null","w",stderr);
#else
    freopen("paint.in","r",stdin);
    freopen("paint.out","w",stdout);
#endif
    ios::sync_with_stdio(false),cin.tie(nullptr),cout.tie(nullptr);
    cin>>n>>(tmp+1);
    RNG(i,1,n) A[i]=tmp[i]-'0';
    k=int(count(A+1,A+n+1,1));
    RNG(u,2,n) cin>>fa[u],G[fa[u]].push_back(u);
    if(n==2){ cout<<500000004<<"\n"; return 0; }
    F[1][0]={1,0,0},F[1][1]={0,1,0};
    RNG(i,1,n-2){
        F[i+1][1]=(F[i][1]*n-F[i-1][1]*(i-1)-X{0,0,1}-F[i-1][0])*qpow(n-i,MOD-2);
        F[i+1][0]=(F[i][0]*n-X{0,0,1}-F[i+1][1]-F[i-1][0]*i)*qpow(n-i-1,MOD-2);
    }
    long p,q;
    {
        auto x=F[n-1][0]*n-X{0,0,1}-F[n-2][0]*(n-1);
        auto y=F[n-1][1]*n-F[n-2][1]*(n-2)-X{0,0,1}-F[n-2][0];
        if(!y.b) swap(x,y);
        assert(y.b);
        auto a=(x.a+MOD-(x.b*y.a%MOD)*qpow(y.b,MOD-2)%MOD)%MOD,b=MOD-x.c+(x.b*y.c%MOD*qpow(y.b,MOD-2)%MOD);
        assert(b);
        p=b*qpow(a,MOD-2)%MOD;
        q=(MOD-y.c+MOD-y.a*p%MOD)*qpow(y.b,MOD-2)%MOD;
    }
    RNG(i,1,n-1){
        RNG(t,0,1) H[i][t]=(p*F[i][t].a+q*F[i][t].b+F[i][t].c)%MOD;
    }
    dfs1(1),dfs2(1,0);
    long ans=0;
    RNG(u,1,n) addto(ans,D[u]*H[k][A[u]]);
    cout<<ans<<"\n";
}