#include <vector>
#include <stack>
#include <queue>
#include <list>
#include <bitset>
#include <set>
#include <map>
#include <unordered_set>
#include <unordered_map>
#include <algorithm>
#include <numeric>
#include <iostream>
#include <iomanip>
#include <string>
#include <chrono>
#include <random>
#include <cmath>
#include <cassert>
#include <climits>
#include <cstring>
#include <cstdlib>
#include <functional>
#include <sstream>

using namespace std;


    const static long long MOD = 1000000007LL;

    vector<long long> pow2;
    vector<long long> F;
    vector<long long> IF;

    long long norm(long long a) {
        a %= MOD;
        if (a < 0) {
            a += MOD;
        }
        return a;
    }

    long long inverse(long long a) {
        long long b = MOD;
        long long u = 0, v = 1;
        a = norm(a);
        while (a != 0) {
            long long d = b / a;
            b -= d * a;
            u -= d * v;
            swap(a, b);
            swap(u, v);
        }
        return norm(u);
    }

    long long powM(long long x, long long e) {
        long long res = 1;
        while (e > 0) {
            if (e & 1) {
                res = res * x % MOD;
            }
            x = x * x % MOD;
            e >>= 1;
        }
        return res;
    }

    void init(int n) {
        F.clear();
        F.resize(n + 1, 1);
        IF = F;
        pow2.clear();
        pow2.resize(n + 1, 1);
        for (int i = 1; i <= n; ++i) {
            pow2[i] = pow2[i - 1] * 2 % MOD;
        }

        for (int i = 1; i <= n; ++i) {
            F[i] = F[i - 1] * i % MOD;
        }
        IF[n] = inverse(F[n]);
        for (int i = n - 1; i > 0; --i) {
            IF[i] = IF[i + 1] * (i + 1) % MOD;
        }
    }

    long long nCr(int n, int r) {
        if (r < 0 || r > n) {
            return 0;
        }
        long long res = F[n] * IF[r] % MOD * IF[n - r] % MOD;
        return res;
    }

int main(int argc, char** argv) {
    ios::sync_with_stdio(false);
    cin.tie(0);


    int n;
    cin >> n;
    int rt = 0;

    vector<vector<int>> graph(n + 1);
    vector<int> pars(n + 1, -1);
    for (int i = 1; i <= n; ++i) {
        cin >> pars[i];
        graph[pars[i]].push_back(i);
        graph[i].push_back(pars[i]);
    }

    vector<int> depths(n + 1, 0);
    vector<int> ord;
    ord.reserve(n + 1);
    ord.push_back(rt);
    for (int i = 0; i < ord.size(); ++i) {
        int u = ord[i];
        for (auto v : graph[u]) {
            if (v == pars[u]) {
                continue;
            }
            depths[v] = depths[u] + 1;
            ord.push_back(v);
        }
    }
    assert(ord.size() == n + 1);
    reverse(ord.begin(), ord.end());

    init(n + 10);

    long long res = 0;

    vector<vector<vector<long long>>> dp(n + 1);

    vector<int> maxDepths(n + 1, 0);
    vector<int> sizes(n + 1, 0);
    
    vector<long long> tmp(3, 0);

    long long inv2 = (MOD + 1) / 2;

    for (auto u : ord) {
        int w = -1;
        maxDepths[u] = depths[u];
        sizes[u] = 1;
        for (auto v : graph[u]) {
            if (v == pars[u]) {
                continue;
            }
            if (w == -1 || maxDepths[v] > maxDepths[w]) {
                w = v;
            }
            maxDepths[u] = max(maxDepths[u], maxDepths[v]);
            sizes[u] += sizes[v];
        }

        if (w >= 0) {
            swap(dp[w], dp[u]);
        }

        dp[u].emplace_back(vector<long long>(3, 0));
        dp[u].back()[0] = dp[u].back()[1] = inv2;

        {
            int d = 0;
            for (auto v : graph[u]) {
                if (v == w || v == pars[u]) {
                    continue;
                }
                d = max(d, (int)dp[v].size());
            }

            for (int i = (int) dp[u].size() - 1 - d; i + 1 < dp[u].size(); ++i) {
                (dp[u][i][0] += dp[u][i][2]) %= MOD;
                dp[u][i][2] = 0;
            }
        }
        for (auto v : graph[u]) {
            if (v == w || v == pars[u]) {
                continue;
            }

            auto& A = dp[u];
            auto& B = dp[v];
            int n = A.size();
            int m = B.size();
            // cout << u << " " << v << " " << n << " " << m << endl;
            for (int i = 0, j = n - (m + 1); i < m; ++j, ++i) {
                fill(tmp.begin(), tmp.end(), 0);
                for (int a = 2; a >= 0; --a) {
                    for (int b = 2; b >= 0; --b) {
                        int na = min(a + (b == 1), 2);
                        (tmp[na] += A[j][a] * B[i][b] % MOD) %= MOD;
                    }
                }
                A[j] = tmp;
            }
        }

        // for (int i = 0; i < dp[u].size(); ++i) {
        //     long long sum = 0;
        //     for (int j = 0; j <= 2; ++j) {
        //         sum += dp[u][i][j];
        //     }
        //     if (sum % MOD != 1) {
        //         cout << "Bad " << u << endl;
        //     }
        // }

        // for (auto& V : dp[u]) {
        //     (V[0] += V[2]) %= MOD;
        //     V[2] = 0;
        // }

        // for (int i = 0; i < dp[u].size(); ++i) {
        //     cout << u << "::: " << i << " ::: " << dp[u][i][1] << " " << dp[u][i][0] <<  endl;
        //     if ((dp[u][i][0] + dp[u][i][1]) % MOD != 1) {
        //         cout << u << " " << depths[u] << " " << pow2[sizes[u]] << " " << graph[u].size() << endl;
        //     }
        //     assert((dp[u][i][0] + dp[u][i][1]) % MOD == 1);
        // }
    } 

    for (int i = 0; i < dp[rt].size(); ++i) {
        // cout << i << " " << dp[rt][i][1] << endl;
        (res += dp[rt][i][1]) %= MOD;
    }
    (res *= pow2[n + 1]) %= MOD;
    cout << res << '\n';

    return 0;
}

Measured memory_bytes: 40632