#include <bits/stdc++.h>
#include <unistd.h>
#define ll long long
#define inf 1000000007
#define inf16 0x3f3f3f3f
#define INF 1000000000000000007LL
#define VI vector<int>
#define VPII vector<pair<int, int> >
#define VLL vector<ll>
#define PII pair<int, int>
#define st first
#define nd second
#define pb push_back
#define mp make_pair
#define eb emplace_back
#define endl '\n'
#define ALL(c) (c).begin(), (c).end()
using namespace std;

#define int ll
const int mod = 1e9+7;
const int N = 2007;
const int inv2 = (mod+1)/2;

int n;

VI G[N];

int sub[N];
int dp[2][N][N][3][2]; //[korzen][ile skladowych][stopien korzenia][parzystosc liczby krawedzi]
int c[N];

void dfs(int v)
{
	sub[v] = 1;

	dp[c[v]][v][1][0][0] = 1;

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

		for(int degv = 2; degv >= 0; --degv)
		{
			for(int sv = sub[v]; sv >= 1; --sv)
			{
				for(int pv = 0; pv < 2; ++pv)
				{
					for(int degu = 2; degu >= 0; --degu)
					{
						for(int su = sub[u]; su >= 1; --su)
						{
							for(int pu = 0; pu < 2; ++pu)
							{
								if(degv<=1 && degu<=1) //uzywam krawedzi
								{
									int ndeg = degv+1;
									int np = (pv+pu+1)%2;
									int ns = sv+su-1;
									int mult = 1;

									if(degv==0 && degu==0)
										mult *= 2;
									else if(degv==1 && degu==1)
										mult *= inv2;

									dp[c[v]^1][v][ns][ndeg][np] += mult*dp[c[v]][v][sv][degv][pv]%mod*dp[c[u]][u][su][degu][pu]%mod;
									dp[c[v]^1][v][ns][ndeg][np] %= mod;
								}

								//nie uzywam krawedzi

								int ndeg = degv;
								int np = (pv+pu)%2;
								int ns = sv+su;

								dp[c[v]^1][v][ns][ndeg][np] += dp[c[v]][v][sv][degv][pv]*dp[c[u]][u][su][degu][pu]%mod;
								dp[c[v]^1][v][ns][ndeg][np] %= mod;
							}
						}
					}

					dp[c[v]][v][sv][degv][pv] = 0;
				}	
			}
		}

		c[v] ^= 1;
		sub[v] += sub[u];
	}
}

int32_t main()
{
	ios_base::sync_with_stdio(0);
	cin.tie(NULL);

	cin >> n;

	for(int i = 2; i <= n; ++i)
	{
		int x;
		cin >> x;
		G[x].pb(i);
	}

	dfs(1);

	int fact = 1;
	int ans = 0;

	for(int i = 1; i <= n; ++i)
	{
		fact = fact*i%mod;

		for(int deg = 0; deg < 3; ++deg)
		{
			ans -= fact*dp[c[1]][1][i][deg][1]%mod;
			ans += fact*dp[c[1]][1][i][deg][0]%mod;
			ans %= mod;
		}
	}

	ans += mod;
	ans %= mod;

	cout << ans << endl;
}