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

#define rep(i, N) for (int i = 0; i < N; i++)
#define pb push_back

typedef long long ll;
typedef pair<int, int> i_i;
const ll INF = LLONG_MAX / 2;
const int MOD = 1e9 + 7;

int N;
ll sz[2010], dp[2010][2010][2];
vector<int> G[2010];

void dfs_sz(int u) {
	sz[u] = 1;
	for (int v: G[u]) dfs_sz(v), sz[u] += sz[v];
}

void dfs(int u) {
	int n = 1;
	vector<vector<int> > a(n + 1, vector<int>(3));
	a[0][0] = 1;
	for (int v: G[u]) {
		dfs(v);
		vector<vector<int> > _a(n + sz[v] + 1, vector<int>(3));
		rep(x, n + 1) rep(y, 3) if (a[x][y])
			rep(dx, sz[v] + 1) rep(dy, 2) if (dp[v][dx][dy]) {
				int _x = x + dx, _y = y + dy;
				if (_x > n + sz[v]) continue;
				if (_y > 2) continue;
				_a[_x][_y] = (_a[_x][_y] + (ll)a[x][y] * dp[v][dx][dy]) % MOD;
			}
		a = _a;
		n += sz[v];
	}
	rep(x, n + 1) {
		dp[u][x][0] = (dp[u][x][0] + a[x][0]) % MOD;
		dp[u][x + 1][1] = (dp[u][x + 1][1] - a[x][0]) % MOD;
		dp[u][x][0] = (dp[u][x][0] + a[x][1] * 2) % MOD;
		dp[u][x + 1][1] = (dp[u][x + 1][1] - a[x][1]) % MOD;
		dp[u][x][0] = (dp[u][x][0] + a[x][2] * 2) % MOD;
	}
}

ll F[2010];

int main() {
	F[0] = 1;
	rep(x, 2001) F[x + 1] = F[x] * (x + 1) % MOD;
	cin >> N;
	for (int u = 1; u < N; u++) {
		int p; cin >> p, p--;
		G[p].pb(u);
	}
	dfs_sz(0), dfs(0);
	/*
	rep(u, N) {
		rep(x, N + 1) cout << u << ' ' << x << ' ' << dp[u][x][0] << endl;
	}
	rep(x, N + 1) cout << x << ' ' << dp[0][x][0] << endl;
	*/
	int ans = 0;
	rep(x, N + 1) {
		ans = (ans + (ll)dp[0][x][0] * F[N - x]) % MOD;
	}
	cout << (ans + MOD) % MOD << endl;
}
