#include <bits/stdc++.h>

using namespace std;

template <class T> void mini(T &l, T r) {l = min(l, r);}
template <class T> void maxi(T &l, T r) {l = max(l, r);}

string to_string(string s) {
  return '"' + s + '"';
}

string to_string(bool a) {
  return a ? "T" : "F";
}
 
string to_string(const char* s) {
  return to_string((string) s);
}
 
template <typename A, typename B>
string to_string(pair<A, B> p) {
  return "(" + to_string(p.first) + ", " + to_string(p.second) + ")";
}
 
template <typename A>
string to_string(A v) {
  bool first = true;
  string res = "{";
  for (const auto &x : v) {
    if (!first) {
      res += ", ";
    }
    first = false;
    res += to_string(x);
  }
  res += "}";
  return res;
}
 
void debug_out() { cerr << endl; }
 
template <typename Head, typename... Tail>
void debug_out(Head H, Tail... T) {
  cerr << " " << to_string(H);
  debug_out(T...);
}
 
#ifdef LOCAL
#define debug(...) cerr << "[" << #__VA_ARGS__ << "]:", debug_out(__VA_ARGS__)
#else
#define debug(...) 42
#endif

#define mid ((l + r) >> 1)

const int N = 2e5 + 10, B = 1e9 + 7;

using ll = long long;

ll mypw(ll a, int b) {
  ll ret = 1;
  while (b) {
    if (b & 1) ret = ret * a % B;
    b >>= 1;
    a = a * a % B;
  }
  return ret;
}

struct Modint {
  ll x;

  Modint(ll a = 0) {
    x = a % B;
    if (x < 0) x += B;
  }

  Modint inv() {
    return Modint(mypw(x, B - 2));
  }

};

Modint operator+(Modint a, Modint b) {
  return Modint(a.x + b.x);
}

Modint operator-(Modint a, Modint b) {
  return Modint(a.x - b.x);
}

Modint operator*(Modint a, Modint b) {
  return Modint(a.x * b.x);
}

Modint operator/(Modint a, Modint b) {
  return a * b.inv();
}

string to_string(Modint a) {
  return to_string(a.x);
}

int n;
int p[N];
Modint ans;
vector <int> e[N];
int sz[N];
int dep[N];
int rt[N];

struct Node {
  int sz;
  Modint p1, p2;
  int ls, rs;
  bool p2_0;
  Node() {
    sz = -1;
    p1 = 0, p2 = 0;
    ls = rs = 0;
    p2_0 = 0;
  }
} pool[N * 40];

struct Seg {

  int node_cnt;

  void init() {
    node_cnt = -1;
    get_node();
  }

  int get_node() {
    pool[++node_cnt] = Node();
    return node_cnt;
  }

  void push_down(int p, int l, int r) {
    if (pool[p].p2_0) {
      int t = pool[p].ls;
      if (t) {
        pool[t].p2_0 = 1;
        pool[t].p2 = 0;
      }
      t = pool[p].rs;
      if (t) {
        pool[t].p2_0 = 1;
        pool[t].p2 = 0;
      }
      pool[p].p2_0 = 0;
    }
  }

  Node query(int &p, int l, int r, int x) {
    if (!p) p = get_node();
    if (l == r) {
      return pool[p];
    }
    push_down(p, l, r);
    if (x <= mid) return query(pool[p].ls, l, mid, x);
    else return query(pool[p].rs, mid + 1, r, x);
  }

  void update(int &p, int l, int r, int x, Modint p1, Modint p2) {
    if (!p) p = get_node();
    pool[p].sz = max(pool[p].sz, x);
    if (l == r) {
      pool[p].p1 = p1;
      pool[p].p2 = p2;
      pool[p].p2_0 = 0;
      return;
    }
    push_down(p, l, r);
    if (x <= mid) update(pool[p].ls, l, mid, x, p1, p2);
    else update(pool[p].rs, mid + 1, r, x, p1, p2);
  }

  void set_p2_0(int &p, int l, int r) {
    if (!p) return;
    pool[p].p2_0 = 1;
  }

  vector < pair <int, Node> > dfs(int &p, int l, int r) {
    vector < pair <int, Node> > ret;
    dfs(p, l, r, ret);
    return ret;
  }

  void dfs(int &p, int l, int r, vector < pair <int, Node> > &v) {
    if (!p) return;
    if (l == r) {
      v.emplace_back(l, pool[p]);
      return;
    }
    push_down(p, l, r);
    dfs(pool[p].ls, l, mid, v);
    dfs(pool[p].rs, mid + 1, r, v);
  }


} seg;

void dfs(int u, int d) {
  sz[u] = 1;
  dep[u] = d;
  for (auto v : e[u]) {
    dfs(v, d + 1);
    sz[u] += sz[v];
  }
}

void mrg(int &a, int &b) {
  if (pool[b].sz > pool[a].sz) swap(a, b);
  vector < pair <int, Node> > v = seg.dfs(b, 0, n);
  for (auto kv : v) {
    //if (kv.second.p2_0) kv.second.y2 = 0;
    int x = kv.first;
    auto y1 = kv.second.p1, y2 = kv.second.p2;
    auto av = seg.query(a, 0, n, x);
    //if (av.p2_0) av.p2 = 0;
    auto a1 = av.p1, a2 = av.p2;
    auto t1 = a1 * (1 - y1 - y2) + (1 - a1 - a2) * y1;
    auto t2 = a1 * (y1 + y2) + a2 + (1 - a1 - a2) * y2;
    seg.update(a, 0, n, x, t1, t2);
  }
}

int dfs(int u) {
  int &p = rt[u];
  seg.update(p, 0, n, dep[u], Modint(1) / 2, 0);
  for (auto v : e[u]) {
    auto t = dfs(v);
    mrg(p, t);
  }
  seg.set_p2_0(p, 0, n);
  return p;
}

void run() {
  cin >> n;
  for (int i = 1; i <= n; ++i) {
    cin >> p[i];
    e[p[i]].push_back(i);
  }
  dfs(0, 0);
  auto t = dfs(0);

  vector < pair <int, Node> > v = seg.dfs(t, 0, n);

  for (auto kv : v) {
    ans = ans + kv.second.p1 * mypw(2, n + 1);
  }

  cout << ans.x << endl;

}

int main() {
  run();
  return 0;
}


Measured memory_bytes: 352748