#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;

long long mod = 1000000007;
long long N, B[1 << 18];
long long sz[2009];
long long dp[2009][4009][2];
long long dr[2009][4009][3];
long long fact[2009];
long long power[2009];
vector<int> X[1 << 18];

long long Div(long long a) {
	if (a % 2 == 1) a += mod;
	return a / 2;
}

void dfs(int pos) {
	sz[pos] = 1;
	for (int i : X[pos]) {
		dfs(i);
		sz[pos] += sz[i];
	}

	// DP の計算
	for (int i = 0; i <= X[pos].size(); i++) {
		for (int j = 0; j <= sz[pos]; j++) {
			dr[i][j][0] = 0; dr[i][j][1] = 0; dr[i][j][2] = 0;
		}
	}

	dr[0][0][0] = 1; int S = 0;
	for (int i = 0; i < X[pos].size(); i++) {
		int to = X[pos][i];
		for (int j = 0; j <= S; j++) {
			for (int k = 0; k <= sz[to]; k++) {
				dr[i + 1][j + k][0] += dr[i][j][0] * dp[to][k][0];
				dr[i + 1][j + k][1] += dr[i][j][1] * dp[to][k][0];
				dr[i + 1][j + k][2] += dr[i][j][2] * dp[to][k][0];
				dr[i + 1][j + k][1] += dr[i][j][0] * dp[to][k][1];
				dr[i + 1][j + k][2] += dr[i][j][1] * dp[to][k][1];
				dr[i + 1][j + k][0] %= mod;
				dr[i + 1][j + k][1] %= mod;
				dr[i + 1][j + k][2] %= mod;
			}
		}
		S += sz[to];
	}

	for (int i = 0; i <= S; i++) {
		dp[pos][i + X[pos].size() - 0][0] += Div(dr[X[pos].size()][i][0]);
		dp[pos][i + X[pos].size() - 0][1] += dr[X[pos].size()][i][0];
		if (X[pos].size() >= 1) {
			dp[pos][i + X[pos].size() - 1][0] -= dr[X[pos].size()][i][1];
			dp[pos][i + X[pos].size() - 1][1] -= dr[X[pos].size()][i][1];
		}
		if (X[pos].size() >= 2) {
			dp[pos][i + X[pos].size() - 2][0] += dr[X[pos].size()][i][2];
		}
	}
	for (int i = 0; i <= N; i++) {
		dp[pos][i][0] = (dp[pos][i][0] + mod * mod) % mod;
		dp[pos][i][1] = (dp[pos][i][1] + mod * mod) % mod;
		dp[pos][i][2] = (dp[pos][i][2] + mod * mod) % mod;
	}
}

void init() {
	power[0] = 1; fact[0] = 1;
	for (int i = 1; i <= 2000; i++) power[i] = (2LL * power[i - 1]) % mod;
	for (int i = 1; i <= 2000; i++) fact[i] = (1LL * i * fact[i - 1]) % mod;
}

int main() {
	cin >> N;
	for (int i = 2; i <= N; i++) cin >> B[i];
	for (int i = 2; i <= N; i++) X[B[i]].push_back(i);

	init();
	dfs(1);
	
	long long Answer = 0;
	for (int i = 1; i <= N; i++) {
		long long v1 = dp[1][i - 1][0];
		long long v2 = fact[i];
		long long v3 = power[i];
		Answer += (v1 * v2 % mod) * v3 % mod;
		Answer %= mod;
	}
	cout << Answer << endl;
	return 0;
}