#include <bits/stdc++.h>
#define syosu(x) fixed<<setprecision(x)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> P;
typedef pair<double,double> pdd;
typedef pair<ll,ll> pll;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef vector<double> vd;
typedef vector<vd> vvd;
typedef vector<ll> vl;
typedef vector<vl> vvl;
typedef vector<string> vs;
typedef vector<P> vp;
typedef vector<vp> vvp;
typedef vector<pll> vpll;
typedef pair<int,P> pip;
typedef vector<pip> vip;
const int inf=1<<30;
const ll INF=1ll<<60;
const double pi=acos(-1);
const double eps=1e-10;
const ll mod=1e9+7;
const int dx[4]={0,1,0,-1},dy[4]={1,0,-1,0};

int n;
vvi g;
vl F;
vvl DP,dp[3];

int dfs(int v){
	int S=1;
	dp[0][v][0]++;
	vp b;
	for(int i=0;i<g[v].size();i++){
		int u=g[v][i];
		b.push_back({dfs(u),u});
	}
	sort(b.rbegin(),b.rend());
	for(int i=0;i<b.size();i++){
		int T=b[i].first,u=b[i].second;
		S+=T;
		vl a[3];
		a[0]=a[1]=a[2]=vl(n);
		a[0][0]++;
		if(!i){
			a[0]=DP[u];
			for(int j=1;j<S;j++) a[1][j]=(dp[0][u][j-1]*2+dp[1][u][j-1])%mod;
		}
		else{
			for(int j=0;j<3;j++) for(int k=1;k<S;k++) for(int l=max(0,k-T);l<=k;l++){
				(a[j][k]+=dp[j][v][l]*DP[u][k-l])%=mod;
				if(j&&l<k) (a[j][k]+=dp[j-1][v][l]*(dp[0][u][k-l-1]*(3-j)+dp[1][u][k-l-1]*(j==1?1:mod/2+1)%mod))%=mod;
			}
		}
		dp[0][v]=a[0];
		dp[1][v]=a[1];
		dp[2][v]=a[2];
	}
	for(int i=0;i<3;i++) for(int j=0;j<S;j++) (DP[v][j]+=dp[i][v][j])%=mod;
	return S;
}

int main(){
	cin>>n;
	g=vvi(n);
	F=vl(n+1,1);
	DP=dp[0]=dp[1]=dp[2]=vvl(n,vl(n));
	for(int i=1;i<n;i++){
		int v;
		cin>>v;
		g[v-1].push_back(i);
	}
	dfs(0);
	ll res=0;
	for(int i=1;i<=n;i++) F[i]=F[i-1]*i%mod;
	for(int i=0;i<n;i++){
		ll tmp=F[n-i]*DP[0][i]%mod;
		if(i%2==1) tmp=(mod-tmp)%mod;
		(res+=tmp)%=mod;
	}
	cout<<res<<endl;
}