#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
typedef long long ll ;
#define rep(i, a, b) for (int i = a; i <= b; ++ i) 
const int N = 4e5 + 5, mo = 1e9 + 7 ;
using namespace std ;

int st[N], top, dp[N][2], f[N], b[N], up[N] ;
int n, e, ter[N], nxt[N], lnk[N], fac[N], cnt[N], dep[N], fa[N][21], a[N], tim, L[N], R[N] ;

void add(int x, int y) {
    ter[++ e] = y, nxt[e] = lnk[x], lnk[x] = e ;
}

void upd(int &x, int y) {
    x = (x + y) % mo ;
}

int mul(int a, int b) {
    int ans = 1, w = a % mo ;
    for ( ; b ; b /= 2) {
        if (b & 1) ans = (ll) ans * w % mo ;
        w = (ll) w * w % mo ;
    }
    return ans ;
}

void dfs(int p) {
    L[p] = ++ tim ;
    ++ cnt[dep[p]] ;
    for (int i = 1 ; (1 << i) <= dep[p] ; ++ i) 
        fa[p][i] = fa[fa[p][i - 1]][i - 1] ;
    for (int i = lnk[p] ; i; i = nxt[i]) {
        fa[ter[i]][0] = p ;
        dep[ter[i]] = dep[p] + 1 ;
        dfs(ter[i]) ;
    }
    R[p] = tim ;
}

int lca(int x, int y) {
    if (dep[x] < dep[y]) swap(x, y) ;
    int res = dep[x] - dep[y] ;
    for (int i = 20 ; i >= 0 ; -- i) if ((res >> i) & 1) 
        x = fa[x][i] ;
    if (x == y) return x ;
    for (int i = 20 ; i >= 0 ; -- i) if (fa[x][i] != fa[y][i])
        x = fa[x][i], y = fa[y][i] ;
    return fa[x][0] ;
}

bool _cmp(int x, int y) {
    return L[x] < L[y] ;
}

int solve(int d) {
    a[0] = 0 ;
    for ( ; top <= n && dep[st[top]] == d ; ++ top) a[++ a[0]] = st[top] ;
    sort(a + 1, a + a[0] + 1, _cmp) ;
    int tmp = a[0] ;
    rep(i, 1, tmp - 1) a[++ a[0]] = lca(a[i], a[i + 1]) ;
    sort(a + 1, a + a[0] + 1, _cmp) ;
    a[0] = unique(a + 1, a + a[0] + 1) - (a + 1) ;
    tmp = 1, b[1] = a[1], up[a[1]] = 0 ;
    rep(i, 2, a[0]) {
        for ( ; !(L[a[i]] >= L[b[tmp]] && L[a[i]] <= R[b[tmp]]); -- tmp) ;
        up[a[i]] = b[tmp] ;
        b[++ tmp] = a[i] ;
    }
    rep(i, 1, a[0]) f[a[i]] = 1, dp[a[i]][1] = 0, dp[a[i]][0] = 1 ;
    for (int i = a[0] ; i; -- i) {
        int u = a[i], v = up[u] ;
        dp[u][1] = (ll) dp[u][1] * f[u] % mo ;
        upd(dp[u][0], mo - dp[u][1]) ;
        if (dep[u] == d) dp[u][0] = dp[u][1] = 1 ;
        if (!v) continue ;
        f[v] = (ll) f[v] * dp[u][0] % mo ;
        upd(dp[v][1], (ll) dp[u][1] * mul(dp[u][0], mo - 2) % mo) ;
        dp[v][0] = (ll) dp[v][0] * (dp[u][0] + dp[u][1]) % mo ;
    }
    return dp[a[1]][1] ;
}

bool cmp(int x, int y) {
    return dep[x] < dep[y] ;
}

int main() {
    scanf("%d", &n) ;
    ++ n ;
    int x, ans = 0 ;
    rep(i, 2, n) {
        scanf("%d", &x) ;
        add(x + 1, i) ;
    }
    fac[0] = 1 ;
    rep(i, 1, n) fac[i] = (ll) fac[i - 1] * 2 % mo ;
    dfs(1) ;
    rep(i, 1, n) st[i] = i ;
    top = 1 ;
    sort(st + 1, st + n + 1, cmp) ;
    rep(i, 0, n - 1) if (cnt[i]) {
        upd(ans, (ll) solve(i) * fac[n - cnt[i]] % mo) ;
    } else break ;
    printf("%d\n", ans) ;
    return 0 ; 
}

/*#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
typedef long long ll ;
#define rep(i, a, b) for (int i = a; i <= b; ++ i) 
const int N = 2005, mo = 1e9 + 7 ;
using namespace std ;

int n, e, ter[N], nxt[N], lnk[N], fac[N], D, cnt, a[N], b[N], c[N], dp[N][2] ;

void add(int x, int y) {
    ter[++ e] = y, nxt[e] = lnk[x], lnk[x] = e ;
}

void upd(int &x, int y) {
    x = (x + y) % mo ;
}

void dfs(int p, int dep) {
    if (dep == D) {
        ++ cnt ;
        dp[p][0] = dp[p][1] = 1 ;
        return ;
    }
    for (int i = lnk[p] ; i; i = nxt[i]) dfs(ter[i], dep + 1) ;
    dp[p][1] = 0 ;
    a[0] = 0 ;
    for (int i = lnk[p] ; i; i = nxt[i]) a[++ a[0]] = ter[i] ;
    b[0] = c[a[0] + 1] = 1 ;
    rep(i, 1, a[0]) b[i] = (ll) b[i - 1] * dp[a[i]][0] % mo ;
    for (int i = a[0] ; i; -- i) 
        c[i] = (ll) c[i + 1] * dp[a[i]][0] % mo ;
    rep(i, 1, a[0]) upd(dp[p][1], (ll) b[i - 1] * c[i + 1] % mo * dp[a[i]][1] % mo) ;
    dp[p][0] = 1 ;
    rep(i, 1, a[0]) dp[p][0] = (ll) dp[p][0] * (dp[a[i]][0] + dp[a[i]][1]) % mo ;
    upd(dp[p][0], mo - dp[p][1]) ;
} 

int main() {
    scanf("%d", &n) ;
    ++ n ;
    int x, ans = 0 ;
    rep(i, 2, n) {
        scanf("%d", &x) ;
        add(x + 1, i) ;
    }
    fac[0] = 1 ;
    rep(i, 1, n) fac[i] = (ll) fac[i - 1] * 2 % mo ;
    rep(i, 1, n) {
        cnt = 0 ;
        D = i ;
        dfs(1, 1) ;
        if (!cnt) break ;
        upd(ans, (ll) dp[1][1] * fac[n - cnt] % mo) ;
    }
    printf("%d\n", ans) ;
    return 0 ; 
}*/

Measured memory_bytes: 46848