#pragma comment(linker, "/STACK:512000000")
#define _CRT_SECURE_NO_WARNINGS
/*#if !defined(YCM) && 0
#define _FORTIFY_SOURCE 0
#pragma GCC optimize("Ofast,no-stack-protector")
#pragma GCC target("avx,tune=native")
#include <stdio.h>
#endif*/
//#include "testlib.h"
#include <bits/stdc++.h>
using namespace std;

#define all(a) a.begin(), a.end()
typedef long long li;
typedef long double ld;
void solve(bool);
void precalc();
clock_t start;
int main() {
#ifdef AIM
    freopen("/home/alexandero/ClionProjects/ACM/input.txt", "r", stdin);
    //freopen("/home/alexandero/ClionProjects/ACM/output.txt", "w", stdout);
    //freopen("out.txt", "w", stdout);
#else
    //freopen("input.txt", "r", stdin);
  //freopen("output.txt", "w", stdout);
#endif
    start = clock();
    int t = 1;
    cout.sync_with_stdio(0);
    cin.tie(0);
    precalc();
    cout.precision(10);
    cout << fixed;
    //cin >> t;
    int testNum = 1;
    while (t--) {
        //cout << "Case #" << testNum++ << ": ";
        //cerr << testNum << endl;
        solve(true);
        //cerr << testNum - 1 << endl;
    }
    cout.flush();
#ifdef AIM1
    while (true) {
      solve(false);
  }
#endif

#ifdef AIM
    cerr << "\n\n time: " << (clock() - start) / 1.0 / CLOCKS_PER_SEC << "\n\n";
#endif

    return 0;
}

//BE CAREFUL: IS INT REALLY INT?

template<typename T>
T binpow(T q, T w, T mod) {
    if (!w)
        return 1 % mod;
    if (w & 1)
        return q * 1LL * binpow(q, w - 1, mod) % mod;
    return binpow(q * 1LL * q % mod, w / 2, mod);
}

template<typename T>
T gcd(T q, T w) {
    while (w) {
        q %= w;
        swap(q, w);
    }
    return q;
}
template<typename T>
T lcm(T q, T w) {
    return q / gcd(q, w) * w;
}

template <typename T>
void make_unique(vector<T>& a) {
    sort(all(a));
    a.erase(unique(all(a)), a.end());
}

template<typename T>
void relax_min(T& cur, T val) {
    cur = min(cur, val);
}

template<typename T>
void relax_max(T& cur, T val) {
    cur = max(cur, val);
}

void precalc() {}

#define int li
const int mod = 1000000007;

const int C = 2010;
int c[C][C];
vector<vector<vector<vector<int>>>> dp;
vector<vector<int>> g;
vector<int> cnt;

void add(int& cur, int val) {
    cur += val;
    cur %= mod;
}

int rev2;

using DP = vector<vector<vector<int>>>;

void dfs(int v) {
    cnt[v] = 1;
    dp[v].assign(1, vector<vector<int>>(2, vector<int>(3, 0)));
    dp[v][0][0][0] = 1;
    DP new_dp;
    int children = 0;
    for (int to : g[v]) {
        dfs(to);
        new_dp.assign(cnt[v] + cnt[to], vector<vector<int>>(2, vector<int>(3, 0)));

        for (int cur_diff = 0; cur_diff < dp[v].size(); ++cur_diff) {
            for (int cur_k = 0; cur_k < 2; ++cur_k) {
                for (int cur_deg = 0; cur_deg <= 2; ++cur_deg) {
                    int cur_dp = dp[v][cur_diff][cur_k][cur_deg];
                    if (!cur_dp) continue;
                    for (int new_diff = 0; new_diff < dp[to].size(); ++new_diff) {
                        for (int new_k = 0; new_k < 2; ++new_k) {
                            for (int new_deg = 0; new_deg <= 2; ++new_deg) {
                                int to_dp = dp[to][new_diff][new_k][new_deg];
                                if (!to_dp) continue;
                                add(new_dp[cur_diff + new_diff][cur_k ^ new_k][cur_deg], cur_dp * to_dp);
                                if (cur_deg == 2 || new_deg == 2) {
                                    continue;
                                }
                                int dp_deg = cur_deg + 1;
                                int dp_k = (cur_k ^ new_k ^ 1);
                                int dp_diff = cur_diff + new_diff + 1;
                                int val = cur_dp * to_dp % mod;
                                if (new_deg && cur_deg) {
                                    val = val * rev2 % mod;
                                } else if (cur_deg == 0 && new_deg == 0) {
                                    val = val * 2 % mod;
                                }
                                add(new_dp[dp_diff][dp_k][dp_deg], val);
                            }
                        }
                    }
                }
            }
        }
        cnt[v] += cnt[to];
        dp[v].swap(new_dp);
        ++children;
    }
}

void solve(bool read) {
    for (int i = 0; i < C; ++i) {
        c[i][i] = c[i][0] = 1;
        for (int j = 1; j < i; ++j) {
            c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mod;
        }
    }
    rev2 = binpow(2LL, mod - 2, mod);
    int n;
    cin >> n;
    g.resize(n);
    cnt.assign(n, 0);
    dp.resize(n);
    for (int i = 1; i < n; ++i) {
        int p;
        cin >> p;
        --p;
        g[p].push_back(i);
    }
    dfs(0);

    vector<int> fact(2 * n + 5, 1);
    for (int i = 1; i < fact.size(); ++i) {
        fact[i] = fact[i - 1] * i % mod;
    }

    int ans = 0;
    for (int diff = 0; diff < n; ++diff) {
        int cur_fact = fact[n - diff];
        for (int par_k = 0; par_k < 2; ++par_k) {

            /*for (int j = 0; j < 3; ++j) {
                if (dp[0][diff][par_k][j] > 0) {
                    cout << diff << " " << par_k << " " << j << " " << dp[0][diff][par_k][j] << endl;
                }
            }*/

            //cout << diff << " " << par_k << " " << cur_fact * (dp[0][diff][par_k][0] + dp[0][diff][par_k][1] + dp[0][diff][par_k][2]) << endl;

            add(ans, cur_fact * (dp[0][diff][par_k][0] + dp[0][diff][par_k][1] + dp[0][diff][par_k][2]));
            cur_fact = -cur_fact;
        }
    }
    if (ans < 0) {
        ans += mod;
    }

    cout << ans << endl;

}