#include<bits/stdc++.h>
using namespace std;

#define int long long

#define rep(i,n) for(int i=0;i<(n);i++)
#define pb push_back
#define all(v) (v).begin(),(v).end()
#define fi first
#define se second
typedef vector<int>vint;
typedef pair<int,int>pint;
typedef vector<pint>vpint;

template<typename A,typename B>inline void chmin(A &a,B b){if(a>b)a=b;}
template<typename A,typename B>inline void chmax(A &a,B b){if(a<b)a=b;}

const int mod=1000000007;
inline void add(int &a,int b){
    a+=b;
    if(a>=mod)a-=mod;
}
int mpow(int n,int m){
    int ret=1;
    while(m){
        if(m&1)ret=ret*n%mod;
        n=n*n%mod;
        m>>=1;
    }
    return ret;
}
const int FACT_SIZE=1111111;
int fact[FACT_SIZE];
int inv[FACT_SIZE];
struct fact_exec{
    fact_exec(){
        fact[0]=1;
        for(int i=1;i<FACT_SIZE;i++)fact[i]=fact[i-1]*i%mod;
        inv[FACT_SIZE-1]=mpow(fact[FACT_SIZE-1],mod-2);
        for(int i=FACT_SIZE-2;i>=0;i--)inv[i]=inv[i+1]*(i+1)%mod;
    }
}factexec;
int nCk(int n,int k){
    return fact[n]*inv[k]%mod*inv[n-k]%mod;
}
int nPk(int n,int k){
    return fact[n]*inv[n-k]%mod;
}

int inv2=(mod+1)/2;

int N;
vint G[2222];

int dp[2222][2222][3];
int sz[2222];
int tmp[2222][3];

void dfs(int v){
    sz[v]=1;
    dp[v][1][0]=1;

    for(auto u:G[v]){
        dfs(u);

        rep(i,sz[v]+sz[u]+1)rep(j,3)tmp[i][j]=0;

        for(int i=1;i<=sz[v];i++){
            for(int j=1;j<=sz[u];j++){
                //0 0
                add(tmp[i+j][0],dp[v][i][0]*dp[u][j][0]%mod);
                add(tmp[i+j-1][1],dp[v][i][0]*dp[u][j][0]*2%mod);

                //0 1
                add(tmp[i+j][0],dp[v][i][0]*dp[u][j][1]%mod);
                add(tmp[i+j-1][1],dp[v][i][0]*dp[u][j][1]%mod);

                //0 2
                add(tmp[i+j][0],dp[v][i][0]*dp[u][j][2]%mod);

                //1 0
                add(tmp[i+j][1],dp[v][i][1]*dp[u][j][0]%mod);
                add(tmp[i+j-1][2],dp[v][i][1]*dp[u][j][0]%mod);

                //1 1
                add(tmp[i+j][1],dp[v][i][1]*dp[u][j][1]%mod);
                add(tmp[i+j-1][2],dp[v][i][1]*dp[u][j][1]%mod*inv2%mod);

                //1 2
                add(tmp[i+j][1],dp[v][i][1]*dp[u][j][2]%mod);

                //2 0
                add(tmp[i+j][2],dp[v][i][2]*dp[u][j][0]%mod);

                //2 1
                add(tmp[i+j][2],dp[v][i][2]*dp[u][j][1]%mod);

                //2 2
                add(tmp[i+j][2],dp[v][i][2]*dp[u][j][2]%mod);
            }
        }

        sz[v]+=sz[u];
        rep(i,sz[v]+1)rep(j,3)dp[v][i][j]=tmp[i][j];
    }
}

signed main(){
    cin>>N;
    for(int i=1;i<N;i++){
        int b;cin>>b;
        b--;
        G[b].pb(i);
    }

    dfs(0);
    int ans=0;
    for(int i=1;i<=N;i++){
        int ei=0;
        rep(j,3)add(ei,dp[0][i][j]);
        ei=ei*fact[i]%mod;

        if((N-i)&1){
            ans=(ans-ei+mod)%mod;
        }
        else{
            ans=(ans+ei)%mod;
        }
    }
    cout<<ans<<endl;
    return 0;
}
