#include<stdio.h>
#include<math.h>
#include<algorithm>
#include<queue>
#include<deque>
#include<string>
#include<string.h>
#include<vector>
#include<set>
#include<map>
#include<bitset>
#include<stdlib.h>
#include<cassert>
using namespace std;
const long long mod=1000000007;
const long long inf=mod*mod;
const long long d2=500000004;
const double EPS=1e-6;
const double PI=acos(-1.0);
int ABS(int a){return max(a,-a);}
long long ABS(long long a){return max(a,-a);}
int g[2100];
int sz[2100];
long long dp[2][2100];
long long dp2[3][2100][2100];
long long fact[3100];
void solve(int a){
	int n=g[a].size()
	dp2[0][0][0]=1;
	for(int i=0;i<g[a].size();i++){
		solve(g[a][i]);
		int T=sz[a];
		sz[a]+=sz[g[a][i]];
		for(int j=0;j<=sz[a];j++){
			dp2[0][i+1][j]=dp2[1][i+1][j]=dp2[2][i+1][j]=0;
		}
		for(int j=0;j<=F;j++){
			for(int k=0;k<=sz[g[a][i]];k++){
				dp2[0][i+1][j+k]=(dp2[0][i+1][j+k]+dp2[0][i][j]*dp[0][k])%mod;
				dp2[1][i+1][j+k]=(dp2[0][i+1][j+k]+dp2[0][i][j]*dp[1][k])%mod;
				dp2[1][i+1][j+k]=(dp2[0][i+1][j+k]+dp2[1][i][j]*dp[0][k])%mod;
				dp2[2][i+1][j+k]=(dp2[0][i+1][j+k]+dp2[1][i][j]*dp[1][k])%mod;
				dp2[2][i+1][j+k]=(dp2[0][i+1][j+k]+dp2[2][i][j]*dp[0][k])%mod;
			}
		}
	}
	for(int i=0;i<=sz[g[a][i]];i++){
		dp[0][i]=(dp2[0][n][i]+dp2[1][n][i]*2+dp2[2][n][i])%mod;
		dp[1][i+1]=dp2[1][n][i]%mod;
	}
	sz[a]++;
}
int main(){
	int a;scanf("%d",&a);
	for(int i=0;i<a-1;i++){
		int p;scanf("%d",&p);
		p--;
		g[p].push_back(i+1);
	}
	solve(0);
	long long ret=0;
	fact[0]=1;
	for(int i=1;i<3100;i++)fact[i]=fact[i-1]*i%mod;

	for(int i=0;i<=a;i++){
		if(i%2==0)ret=(ret+dp[0][i]*fact[a-i]%mod)%mod;
		else ret=(ret+mod-dp[0][i]*fact[a-i]%mod)%mod;
	}
	printf("%lld\n",ret);
}