#include <bits/stdc++.h>
using namespace std;
const long long mod = 1000000007;
const long long inf = mod * mod;
const long long d2 = 500000004;
const double EPS = 1e-6;
const double PI = acos(-1.0);
int ABS(int a) { return max(a, -a); }
long long ABS(long long a) { return max(a, -a); }
vector<int> g[2100];
int sz[2100];
long long dp[2][2100][2100];
long long dp2[3][2100][4100];
long long fact[3100];
void solve(int a) {
  int n = g[a].size();
  dp2[0][0][0] = 1;
  for (int i = 0; i < g[a].size(); i++) {
    solve(g[a][i]);
    int T = sz[a];
    sz[a] += sz[g[a][i]];
    int to = g[a][i];
    for (int j = 0; j <= sz[a] + 10; j++) {
      dp2[0][i + 1][j] = dp2[1][i + 1][j] = dp2[2][i + 1][j] = 0;
    }
    for (int j = 0; j <= T; j++) {
      for (int k = 0; k <= sz[g[a][i]] + 1; k++) {
        dp2[0][i + 1][j + k] =
            (dp2[0][i + 1][j + k] + dp2[0][i][j] * dp[0][to][k]) % mod;
        dp2[1][i + 1][j + k] =
            (dp2[1][i + 1][j + k] + dp2[0][i][j] * dp[1][to][k]) % mod;
        dp2[1][i + 1][j + k] =
            (dp2[1][i + 1][j + k] + dp2[1][i][j] * dp[0][to][k]) % mod;
        dp2[2][i + 1][j + k] =
            (dp2[2][i + 1][j + k] + dp2[1][i][j] * dp[1][to][k]) % mod;
        dp2[2][i + 1][j + k] =
            (dp2[2][i + 1][j + k] + dp2[2][i][j] * dp[0][to][k]) % mod;
      }
    }
  }
  for (int i = 0; i <= sz[a]; i++) {
    dp[0][a][i] = (dp2[0][n][i] + dp2[1][n][i] * 2 + dp2[2][n][i] * 2) % mod;
    dp[1][a][i + 1] = (dp2[0][n][i] + dp2[1][n][i]) % mod;
  }
  sz[a]++;
}
int main() {
  int a;
  scanf("%d", &a);
  for (int i = 0; i < a - 1; i++) {
    int p;
    scanf("%d", &p);
    p--;
    g[p].push_back(i + 1);
  }
  solve(0);
  long long ret = 0;
  fact[0] = 1;
  for (int i = 1; i < 3100; i++) fact[i] = fact[i - 1] * i % mod;
  for (int i = 0; i <= a; i++) {
    if (i % 2 == 0)
      ret = (ret + dp[0][0][i] * fact[a - i] % mod) % mod;
    else
      ret = (ret + mod - dp[0][0][i] * fact[a - i] % mod) % mod;
  }
  printf("%lld\n", ret);
}
