#include <bits/stdc++.h>
#define show(x) std::cerr << #x << " = " << (x) << std::endl
using ll = long long;
template <typename T>
constexpr T MOD = static_cast<T>(1000000007);
template <typename T>
constexpr T INF = std::numeric_limits<T>::max() / 16;
template <typename T, T mod = MOD<T>>
class ModInt
{
private:
    T value;
    static T pow(const T p, const T n) { return p < 0 ? pow(mod + p, n) : n == 0 ? 1 : n % 2 == 1 ? pow(p, n - 1) * p % mod : pow(p * p % mod, n / 2); }
    static T inv(const T p) { return pow(p, mod - 2); }

public:
    ModInt() : value{0} {}
    ModInt(const T val) : value{((val % mod) + mod) % mod} {}
    ModInt(const ModInt<T, mod>& n) : value{n()} {}
    ModInt<T, mod>& operator=(const ModInt<T, mod>& n) { return value = n(), (*this); }
    ModInt<T, mod>& operator=(const T v) { return value = (mod + v % mod) % mod, (*this); }
    ModInt<T, mod> operator+() const { return *this; }
    ModInt<T, mod> operator-() const { return ModInt{mod - value}; }
    ModInt<T, mod> operator+(const ModInt<T, mod>& val) const { return ModInt{value + val()}; }
    ModInt<T, mod> operator-(const ModInt<T, mod>& val) const { return ModInt{value - val() + mod}; }
    ModInt<T, mod> operator*(const ModInt<T, mod>& val) const { return ModInt{value * val()}; }
    ModInt<T, mod> operator/(const ModInt<T, mod>& val) const { return ModInt{value * inv(val())}; }
    ModInt<T, mod>& operator+=(const ModInt<T, mod>& val) { return (((value += val()) %= mod) += mod) %= mod, (*this); }
    ModInt<T, mod>& operator-=(const ModInt<T, mod>& val) { return (((value -= val()) %= mod) += mod) %= mod, (*this); }
    ModInt<T, mod>& operator*=(const ModInt<T, mod>& val) { return (((value *= val()) %= mod) += mod) %= mod, (*this); }
    ModInt<T, mod>& operator/=(const ModInt<T, mod>& val) { return (((value *= inv(val())) %= mod) += mod) %= mod, (*this); }
    ModInt<T, mod> operator+(const T val) const { return ModInt{value + val}; }
    ModInt<T, mod> operator-(const T val) const { return ModInt{value - val}; }
    ModInt<T, mod> operator*(const T val) const { return ModInt{value * val}; }
    ModInt<T, mod> operator/(const T val) const { return ModInt{value * inv(val)}; }
    ModInt<T, mod>& operator+=(const T val) { return (((value += val) %= mod) += mod) %= mod, (*this); }
    ModInt<T, mod>& operator-=(const T val) { return (((value -= val) %= mod) += mod) %= mod, (*this); }
    ModInt<T, mod>& operator*=(const T val) { return (((value *= val) %= mod) += mod) %= mod, (*this); }
    ModInt<T, mod>& operator/=(const T val) { return (((value *= inv(val)) %= mod) += mod) %= mod, (*this); }
    bool operator==(const ModInt<T, mod>& val) const { return value == val.value; }
    bool operator!=(const ModInt<T, mod>& val) const { return not(*this == val); }
    bool operator==(const T val) const { return value == (mod + val % mod) % mod; }
    bool operator!=(const T val) const { return not(*this == val); }
    T operator()() const { return value; }
};
template <typename T, T mod>
inline ModInt<T, mod> operator+(const T val, const ModInt<T, mod>& n) { return ModInt<T, mod>{n() + val}; }
template <typename T, T mod>
inline ModInt<T, mod> operator-(const T val, const ModInt<T, mod>& n) { return ModInt<T, mod>{-n() + val}; }
template <typename T, T mod>
inline ModInt<T, mod> operator*(const T val, const ModInt<T, mod>& n) { return ModInt<T, mod>{n() * val}; }
template <typename T, T mod>
inline ModInt<T, mod> operator/(const T val, const ModInt<T, mod>& n) { return ModInt<T, mod>(val) / n; }
template <typename T, T mod>
inline bool operator==(const T val, const ModInt<T, mod>& n) { return n == val; }
template <typename T, T mod>
inline bool operator!=(const T val, const ModInt<T, mod>& n) { return not(val == n); }
template <typename T, T mod>
inline std::istream& operator>>(std::istream& is, ModInt<T, mod>& n)
{
    T v;
    return is >> v, n = v, is;
}
template <typename T, T mod>
std::ostream& operator<<(std::ostream& os, const ModInt<T, mod>& n) { return (os << n()); }
using mint = ModInt<int, MOD<int>>;
using mll = ModInt<ll, MOD<ll>>;
using Poly = std::vector<mll>;
Poly operator*(const Poly& p1, const Poly& p2)
{
    const int d1 = p1.size() - 1, d2 = p2.size() - 1, d = d1 + d2;
    Poly ans(d + 1, 0);
    for (int i = 0; i <= d1; i++) {
        for (int j = 0; j <= d2; j++) { ans[i + j] += p1[i] * p2[j]; }
    }
    return ans;
}
Poly operator*(const ll s, const Poly& p2)
{
    const int d = p2.size() - 1;
    Poly ans(d + 1, 0);
    for (int i = 0; i <= d; i++) { ans[i] = p2[i] * s; }
    return ans;
}
Poly operator+(const Poly& p1, const Poly& p2)
{
    const int d1 = p1.size() - 1, d2 = p2.size() - 1, d = std::max(d1, d2);
    Poly ans(d + 1, 0);
    for (int i = 0; i <= d1; i++) { ans[i] += p1[i]; }
    for (int i = 0; i <= d2; i++) { ans[i] += p2[i]; }
    return ans;
}
template <typename T, typename A>
std::ostream& operator<<(std::ostream& os, const std::vector<T, A>& v)
{
    os << "[";
    for (const auto& e : v) { os << e << ","; }
    return (os << "]" << std::endl);
}
int main()
{
    int N;
    std::cin >> N;
    std::vector<std::vector<int>> g(N);
    for (int i = 1, p; i < N; i++) { std::cin >> p, g[p - 1].push_back(i); }
    // vi以下のパス分割全体で、根viの次数がj, 連結成分数k
    // 2^{連結成分数}/2^{孤立点数} の和
    std::vector<std::vector<Poly>> dp(N, std::vector<std::vector<mll>>(3));
    auto dfs = [&](auto&& self, const int s) -> void {
        for (const int to : g[s]) { self(self, to); }
        Poly p0{1}, p1{0}, p2{0};
        for (const int to : g[s]) {
            const Poly z = dp[to][0] + dp[to][1] + dp[to][2];
            const Poly o = 2 * dp[to][0] + dp[to][1];
            auto p0_ = p0, p1_ = p1, p2_ = p2;
            p0 = p0_ * z, p1 = p1_ * z + p0_ * o, p2 = p2_ * z + p1_ * o;
        }
        // show(s);
        // show(p0), show(p1), show(p2);
        dp[s][0] = p0 * Poly{0, 1};
        dp[s][1] = p1;
        for (int i = 1; i < p2.size(); i++) { p2[i - 1] = p2[i] * (MOD<ll> + 1) / 2; }
        p2[p2.size() - 1] = 0;
        dp[s][2] = p2;
    };
    dfs(dfs, 0);
    //    show(dp);
    mll p = 1, ans = 0;
    const Poly f = dp[0][0] + dp[0][1] + dp[0][2];
    for (int i = 0; i <= N; i++) {
        ans += ((N - i) % 2 == 0 ? p * f[i] : (-p) * f[i]);
        p *= (i + 1);
    }
    std::cout << ans << std::endl;
    return 0;
}

Measured memory_bytes: 65280