#include <bits/stdc++.h>
using namespace std;
const int inf = 1 << 30;
const long long INF = 1ll << 60;
const double pi = acos(-1);
const double eps = 1e-10;
const long long mod = 1e9 + 7;
const int dx[4] = {0, 1, 0, -1}, dy[4] = {1, 0, -1, 0};
int n;
vector<vector<int> > g;
vector<long long> F;
vector<vector<long long> > DP, dp[3];
int dfs(int v) {
  int S = 1;
  dp[0][v][0]++;
  for (int i = 0; i < g[v].size(); i++) {
    int u = g[v][i];
    int T = dfs(u);
    S += T;
    vector<long long> a[3];
    a[0] = a[1] = a[2] = vector<long long>(n);
    a[0][0]++;
    if (!i) {
      a[0] = DP[u];
      for (int j = 1; j < S; j++)
        a[1][j] = (dp[0][u][j - 1] * 2 + dp[1][u][j - 1]) % mod;
    } else {
      for (int j = 0; j < 3; j++)
        for (int k = 1; k < S; k++)
          for (int l = max(0, k - T); l <= k; l++) {
            (a[j][k] += dp[j][v][l] * DP[u][k - l]) %= mod;
            if (j && l < k)
              (a[j][k] +=
               dp[j - 1][v][l] *
               (dp[0][u][k - l - 1] * (3 - j) +
                dp[1][u][k - l - 1] * (j == 1 ? 1 : mod / 2 + 1) % mod)) %= mod;
          }
    }
    dp[0][v] = a[0];
    dp[1][v] = a[1];
    dp[2][v] = a[2];
  }
  for (int i = 0; i < 3; i++)
    for (int j = 0; j < S; j++) (DP[v][j] += dp[i][v][j]) %= mod;
  return S;
}
int main() {
  cin >> n;
  g = vector<vector<int> >(n);
  F = vector<long long>(n + 1, 1);
  DP = dp[0] = dp[1] = dp[2] =
      vector<vector<long long> >(n, vector<long long>(n));
  for (int i = 1; i < n; i++) {
    int v;
    cin >> v;
    g[v - 1].push_back(i);
  }
  dfs(0);
  long long res = 0;
  for (int i = 1; i <= n; i++) F[i] = F[i - 1] * i % mod;
  for (int i = 0; i < n; i++) {
    long long tmp = F[n - i] * DP[0][i] % mod;
    if (i % 2 == 1) tmp = (mod - tmp) % mod;
    (res += tmp) %= mod;
  }
  cout << res << endl;
}
