#include <bits/stdc++.h>
using namespace std;
typedef long long lint;
typedef pair<int, int> pi;
const int MAXN = 2005;
const int mod = 1e9 + 7;

int n, sz[MAXN], p[MAXN];
vector<int> gph[MAXN];

lint dp[MAXN][MAXN][2];
lint fact[MAXN];
lint aux[MAXN][3], nxt[MAXN][3];

int main(){
	fact[0] = 1;
	for(int i=1; i<MAXN; i++) fact[i] = fact[i-1] * i % mod;
	cin >> n;
	for(int i=2; i<=n; i++){
		cin >> p[i];
		gph[p[i]].push_back(i);
	}
	for(int i=n; i; i--){
		memset(aux, 0, sizeof(aux));
		aux[0][0] = 1;
		for(auto &j : gph[i]){
			memset(nxt, 0, sizeof(nxt));
			for(int k=0; k<=sz[j]; k++){
				for(int l=0; l<=sz[i]; l++){
					for(int m=0; m<3; m++){
						nxt[k + l][m] += dp[j][k][0] * aux[l][m] % mod;
						if(m < 2) nxt[k + l + 1][m+1] += dp[j][k][1] * aux[l][m] % mod;
					}
				}
			}
			sz[i] += sz[j];
			for(int j=0; j<=sz[i]; j++){
				for(int k=0; k<3; k++) aux[j][k] = nxt[j][k] % mod;
			}
		}
		for(int j=0; j<=sz[i]; j++){
			dp[i][j][0] = aux[j][0] + aux[j][1] * 2 + aux[j][2] * 2;
			dp[i][j][1] = aux[j][0] + aux[j][1];
			dp[i][j][0] %= mod;
			dp[i][j][1] %= mod;
		}
		sz[i]++;
	}
	lint ret = 0;
	for(int i=0; i<=n; i++){
		lint cur = dp[1][i][0] * fact[n - i] % mod;
		if(i & 1) cur = mod - cur;
		ret += cur;
	}
	cout << ret % mod << endl;
}
