#include <bits/stdc++.h>
using namespace std;
const int MN = 2005;
int mod = 1e9 + 7;
int n, sz[MN], p[MN];
vector<int> gph[MN];
int dp[MN][MN][2];
int f[MN];
int a[MN][3], nxt[MN][3];
int main() {
  f[0] = 1;
  for (int i = 1; i < MN; i++) f[i] = f[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(a, 0, sizeof(a));
    a[0][0] = 1;
    for (auto &j : gph[i]) {
      memset(nxt, 0, sizeof(nxt));
      for (int k = 0; k < sz[j] + 1; k++)
        for (int l = 0; l < sz[i] + 1; l++)
          for (int m = 0; m < 3; m++) {
            nxt[k + l][m] += dp[j][k][0] * a[l][m] % mod;
            if (m < 2) nxt[k + l + 1][m + 1] += dp[j][k][1] * a[l][m] % mod;
          }
      sz[i] += sz[j];
      for (int j = 0; j < sz[i] + 1; j++)
        for (int k = 0; k < 3; k++) {
          a[j][k] = nxt[j][k] % mod;
        }
    }
    for (int j = 0; j < sz[i] + 1; j++) {
      dp[i][j][0] = a[j][0] + a[j][1] * 2 + a[j][2] * 2;
      dp[i][j][1] = a[j][0] + a[j][1];
      dp[i][j][0] %= mod;
      dp[i][j][1] %= mod;
    }
    sz[i]++;
  }
  int ret = 0;
  for (int i = 0; i <= n; i++) {
    int cur = dp[1][i][0] * f[n - i] % mod;
    if (i & 1) cur = mod - cur;
    ret += cur;
  }
  cout << ret % mod << endl;
}
