#include<bits/stdc++.h>
using namespace std;

#define TRACE(x) x
#define WATCH(x) TRACE(cout << #x" = " << x << endl)
#define WATCHR(a, b) TRACE(for (auto it=a; it!=b;) cout << *(it++) << " "; cout << endl)
#define WATCHC(V) TRACE({cout << #V" = "; WATCHR(V.begin(), V.end());})

#define sz(x) int((x).size())
#define all(x) (x).begin(), (x).end()

using ll = long long;
using vi = vector<int>;
using vvi = vector<vi>;
using vll = vector<ll>;
using vvll = vector<vll>;
using vb = vector<bool>;
using vs = vector<string>;

mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());

template<int MOD> struct modnum {
    int v;
    modnum() : v(0) {}
    modnum(ll _v) : v(_v % MOD) { if (v < 0) v += MOD; }
    explicit operator int() const { return v; }
    friend istream& operator >> (istream& i, modnum& n) { ll v; i >> v; n = modnum(v); return i; }
    friend ostream& operator << (ostream& o, const modnum& n) { return o << n.v; }

    friend bool operator == (const modnum& a, const modnum& b) { return a.v == b.v; }
    friend bool operator != (const modnum& a, const modnum& b) { return a.v != b.v; }

    modnum& operator += (const modnum& o) { v += o.v; if (v >= MOD) v -= MOD; return *this; }
    modnum& operator -= (const modnum& o) { v -= o.v; if (v < 0) v += MOD; return *this; }
    modnum& operator *= (const modnum& o) { v = int(ll(v) * ll(o.v) % MOD); return *this; }
    modnum operator - () { modnum res; if (v) res.v = MOD - v; return res; }
    friend modnum operator + (const modnum& a, const modnum& b) { return modnum(a) += b; }
    friend modnum operator - (const modnum& a, const modnum& b) { return modnum(a) -= b; }
    friend modnum operator * (const modnum& a, const modnum& b) { return modnum(a) *= b; }

    modnum pow(ll e) const {
        if (e == 0) return 1;
        if (e & 1) return *this * this->pow(e-1);
        return (*this * *this).pow(e/2);
    }

    modnum inv() const {
        int g = MOD, x = 0, y = 1;
        for (int r = v; r != 0; ) {
            int q = g / r;
            g %= r; swap(g, r);
            x -= q * y; swap(x, y);
        }

        assert(g == 1);
        assert(y == MOD || y == -MOD);
        return x < 0 ? x + MOD : x;
    }
    modnum& operator /= (const modnum& o) { return (*this) *= o.inv(); }
    friend modnum operator / (const modnum& a, const modnum& b) { return modnum(a) /= modnum(b); }

    static int totient() {
        int tot = MOD, tmp = MOD;
        for (int p = 2; p * p <= tmp; p++) if (tmp % p == 0) {
            tot = tot / p * (p - 1);
            while (tmp % p == 0) tmp /= p;
        }
        if (tmp > 1) tot = tot / tmp * (tmp - 1);
        return tot;
    }

    static int primitive_root() {
        if (MOD == 1) return 0;
        if (MOD == 2) return 1;

        int tot = totient(), tmp = tot;
        vi tot_pr;
        for (int p = 2; p * p <= tmp; p++) if (tot % p == 0) {
            tot_pr.push_back(p);
            while (tmp % p == 0) tmp /= p;
        }
        if (tmp > 1) tot_pr.push_back(tmp);

        for (int r = 2; r < MOD; r++) if (__gcd(r, MOD) == 1) {
            bool root = true;
            for (int p : tot_pr) root &= modnum(r).pow(tot / p) != 1;
            if (root) return r;
        }
        assert(false);
    }

    static modnum generator() { static modnum g = primitive_root(); return g; }
    static int discrete_log(modnum v) {
        static const int M = ceil(sqrt(MOD));
        static unordered_map<int, int> table;
        if (table.empty()) {
            modnum e = 1;
            for (int i = 0; i < M; i++) { table[e.v] = i; e *= generator(); }
        }
        static modnum f = generator().pow(totient() - M);

        for (int i = 0; i < M; i++) {
            if (table.count(v.v)) return table[v.v] + i * M;
            v *= f;
        }
        assert(false);
    }

    static modnum fact(int n) {
        static vector<modnum<MOD>> fact = { 1 };
        for (assert(n >= 0); fact.size() <= n; )
            fact.push_back(fact.back() * fact.size());
        return fact[n];
    }

    static modnum finv(int n) {
        static vector<modnum<MOD>> finv = { 1 };
        for (assert(n >= 0); finv.size() <= n; )
            finv.push_back(finv.back() / finv.size());
        return finv[n];
    }

    static modnum ncr(int n, int r) {
        if (r < 0 || n < r) return 0;
        return fact(n) * finv(r) * finv(n - r);
    }
};
using mn = modnum<int(1e9 + 7)>;

/*
 * Supports O(1) range queries on an immutable array of items.
 * Requires that type T's combine operator is associate, commutative,
 * and idempotent.
 */
template<typename T> struct sparse_table {
    int S, L;
    vector<T> table;
    T& entry(int l, int i) { return table[l * S + i]; }

    // Constructs the table in O(S log S) time and memory
    sparse_table(const vector<T>& elts = {}) : S(elts.size()) {
        if (S) L = 32 - __builtin_clz(S);
        table.resize(L * S);
        copy(all(elts), table.begin());

        for (int l = 0; l + 1 < L; l++) {
            for (int i = 0; i < S; i++) {
                entry(l + 1, i) = entry(l, i);
                if (i + (1 << l) < S)
                    entry(l + 1, i) = entry(l + 1, i) + entry(l, i + (1 << l));
            }
        }
    }

    // Returns elts[i] + ... + elts[j] in O(1)
    T query(int i, int j) {
        int l = 31 - __builtin_clz(j - i + 1);
        return entry(l, i) + entry(l, j - (1 << l) + 1);
    }
};

/*
 * Supports O(1) lowest common ancestor queries on an immutable tree.
 */
struct lowest_common_ancestor {
    struct visit {
        int node, depth;
        visit operator + (const visit& o) const {
            return depth <= o.depth ? *this : o;
        }
    };

    vi depth, pos;
    sparse_table<visit> euler;

    // Constructs the lookup table in O(V log V) time and memory
    lowest_common_ancestor(const vvi& tree = {}, int root = 0) {
        depth.resize(tree.size());
        pos.resize(tree.size());
        vector<visit> tour;
        if (tree.size()) dfs(tree, tour, root, root);
        euler = sparse_table<visit>(tour);
    }

    void dfs(const vvi& tree, auto& tour, int loc, int par) {
        pos[loc] = tour.size();
        tour.push_back({ loc, depth[loc] });
        for (int nbr : tree[loc]) if (nbr != par) {
            depth[nbr] = depth[loc] + 1;
            dfs(tree, tour, nbr, loc);
            tour.push_back({ loc, depth[loc] });
        }
    }

    // Returns the lowest common ancestor of u and v in O(1)
    int lca(int u, int v) {
        u = pos[u], v = pos[v];
        if (u > v) swap(u, v);
        return euler.query(u, v).node;
    }

    int dist(int u, int v) {
        return depth[u] + depth[v] - 2 * depth[lca(u, v)];
    }

    bool on_path(int u, int v, int inx) {
        return dist(u, v) == (dist(u, inx) + dist(inx, v));
    }
};

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    int N;
    cin >> N;

    vvi layer(N + 1);
    layer[0] = { 0 } ;

    vi par(N + 1), dep(N + 1);
    vvi tree(N + 1);

    for (int i = 1; i <= N; i++) {
        cin >> par[i];
        tree[par[i]].push_back(i);
        dep[i] = dep[par[i]] + 1;
        layer[dep[i]].push_back(i);
    }

    lowest_common_ancestor lca(tree);

    mn ans = 0;

    vector<mn> zero(N + 1), one(N + 1);
    for (auto& ly : layer) if (sz(ly)) {
        set<pair<int, int>> elts;
        for (int v : ly) {
            elts.insert({ lca.pos[v], v });
            zero[v] = 1;
            one[v] = 1;
        }

        priority_queue<pair<int, int>> merge;
        for (auto it = elts.begin(); next(it) != elts.end(); it++) {
            merge.push({ lca.depth[lca.lca(it->second, next(it)->second)], it->second });
        }
        while (!merge.empty()) {
            int d, v; tie(d, v) = merge.top(); merge.pop();

            auto it = elts.find({ lca.pos[v], v });
            if (it == elts.end()) continue;

            auto it2 = next(it);
            if (it2 == elts.end()) continue;

            int meet = lca.lca(it->second, it2->second);
            if (lca.depth[meet] != d) continue;

            while (it != elts.begin() && lca.lca(prev(it)->second, it->second) == meet)
                --it;
            while (it2 != elts.end() && lca.lca(it->second, it2->second) == meet)
                it2++;

            vector<mn> zeros, ones;
            for (auto tmp = it; tmp != it2; elts.erase(tmp++)) {
                zeros.push_back(zero[tmp->second]);
                ones.push_back(one[tmp->second]);
            }

            int T = sz(zeros);
            vector<mn> zpref(T), zsuff(T);
            zpref[0] = 1;
            for (int i = 0; i + 1 < T; i++)
                zpref[i + 1] = zpref[i] * zeros[i];

            zsuff[T - 1] = 1;
            for (int i = T - 1; i > 0; i--)
                zsuff[i - 1] = zsuff[i] * zeros[i];

            mn tot = 1; one[meet] = 0;
            for (int i = 0; i < T; i++) {
                tot *= zeros[i] + ones[i];
                one[meet] += ones[i] * zpref[i] * zsuff[i];
            }
            zero[meet] = tot - one[meet];

            elts.insert({ lca.pos[meet], meet });
            auto put = elts.find({ lca.pos[meet], meet });
            if (put != elts.begin()) {
                merge.push({ lca.depth[lca.lca(prev(put)->second, meet)], prev(put)->second });
            }
            if (next(put) != elts.end()) {
                merge.push({ lca.depth[lca.lca(meet, next(put)->second)], meet });
            }
        }

        assert(sz(elts) == 1);
        ans += one[get<1>(*elts.begin())] * mn(2).pow(N + 1 - sz(ly));
    }

    cout << ans << endl;

    return 0;
}



Measured memory_bytes: 95212