import std.algorithm;
import std.array;
import std.ascii;
import std.bigint;
import std.complex;
import std.container;
import std.conv;
import std.functional;
import std.math;
import std.range;
import std.stdio;
import std.string;
import std.typecons;

auto readInts() {
	return array(map!(to!int)(readln().strip().split()));
}
auto readInt() {
	return readInts()[0];
}
auto readLongs() {
	return array(map!(to!long)(readln().strip().split()));
}
auto readLong() {
	return readLongs()[0];
}

void readlnTo(T...)(ref T t) {
    auto s = readln().split();
    assert(s.length == t.length);
    foreach(ref ti; t) {
        ti = s[0].to!(typeof(ti));
        s = s[1..$];
    }
}

class UnionFindTree {
public:
    this(size_t n) {
        parent = iota(n).array();
        rank = new size_t[](n);
    }
    size_t find(size_t n) {
        if(n == parent[n]) {
            return n;
        } else {
            parent[n] = find(parent[n]);
            return parent[n];
        }
    }
    void unite(size_t n, size_t m) {
        auto p1 = find(n);
        auto p2 = find(m);
        if(n == m) {
            return;
        }
        if(rank[p1] < rank[p2]) {
            swap(p1, p2);
        }
        if(rank[p1] == rank[p2]) {
            ++rank[p1];
        }
        parent[p2] = p1;
    }
private:
    size_t[] parent;
    size_t[] rank;
}
unittest {
    auto uft = new UnionFindTree(10);
    foreach(i; iota(10)) {
        assert(uft.find(i) == i);
    }
    foreach(i; iota(10).filter!(a => a%3 == 0).drop(1)) {
        uft.unite(i, i-3);
    }
    foreach(i; iota(10).filter!(a => a%3 == 0)) {
        assert(uft.find(0) == uft.find(i));
    }
    foreach(i; iota(10).filter!(a => a%3 != 0)) {
        assert(uft.find(0) != uft.find(i));
    }
}

const real eps = 1e-10;

const long p = 1_000_000_000 + 7;

long pow(long n, long k) {
    if(k == 0) {
        return 1;
    }
    long tmp = pow(n, k/2);
    tmp = (tmp * tmp) % p;
    if(k%2 == 1) {
        tmp *= n;
        tmp %= p;
    }
    return tmp;
}

void main(){
    auto n = readLong();
    auto P = readLongs();
    auto g = new long[][](n+1, 0);
    foreach(i; iota(n)) {
        g[P[i]] ~= i+1;
    }

    auto size = new long[](n+1);
    long dfs_size(long i) {
        size[i] = g[i].map!(dfs_size).sum(0L) + 1;
        return size[i];
    }
    dfs_size(0);

    auto height = new long[](n+1);
    long dfs_height(long i) {
        height[i] = reduce!max(0L, g[i].map!(dfs_height)) + 1;
        return height[i];
    }
    dfs_height(0);

    auto depthCount = new long[](n+2);
    void dfs_depthCount(long i, long depth) {
        ++depthCount[depth];
        foreach(v; g[i]) {
            dfs_depthCount(v, depth+1);
        }
    }
    dfs_depthCount(0, 1);

    auto depthCountSum = new long[](n+2);
    foreach(i; iota(1, n+2)) {
        depthCountSum[i] = depthCountSum[i-1] + depthCount[i];
    }

    long[][long] dfs_solve(long i, long depth) {
        if(g[i].empty()) {
            return [depth: [1L, 1L]];
        }

        if(g[i].length == 1) {
            auto ret = dfs_solve(g[i][0], depth+1);
            ret[depth] = [pow(2, size[i] - 1), pow(2, size[i]-1)];
            return ret;
        }

        auto heighest = reduce!((a, b) => height[a] > height[b] ? a: b)(g[i]);
        auto ret = dfs_solve(heighest, depth + 1);
        auto tmp = g[i].map!(a => dfs_solve(a, depth+1))();
        auto h2 = g[i].map!(a => height[a]).array.sort!("a > b")[1];
        foreach(h; iota(1, h2+1)) {
            long[] ones;
            long[] zeros;
            foreach(tmpi; tmp) {
                if((depth + h) in tmpi) {
                    zeros ~= tmpi[depth + h][0];
                    ones ~= tmpi[depth + h][1];
                }
            }
            long zPI = 1;
            foreach(zi; zeros) {
                zPI *= zi;
                zPI %= p;
            }
            long o;
            foreach(x; zip(zeros, ones)) {
                o += (x[1] * pow(x[0], p-2)) % p;
                o %= p;
            }
            o *= zPI;
            o %= p;
            long all = 1;
            foreach(x; zip(zeros, ones)) {
                all *= (x[0] + x[1]) % p;
                all %= p;
            }
            ret[depth+h] = [(all + p - o) % p, o];
        }
        ret[depth] = [pow(2, size[i] - 1), pow(2, size[i] - 1)];
        return ret;
    }

    auto tmp = dfs_solve(0, 1);
    long ans;
    foreach(depth, x; tmp) {
        ans += (x[1] * pow(2, depthCountSum[depth-1])) % p;
        ans %= p;
    }

    writeln(ans);
}

