#include<cstdio>
#include<vector>
#include<iostream>
using namespace std;
const int N=2e5+5,p=1e9+7;
typedef long long ll;
struct node {
    int a,b,c;
    void operator +=(const node &x){
        c=((ll)(a+b+c)*x.c+(ll)(b+c)*x.b+c*(ll)x.a)%p;
        b=((ll)a*x.b+(ll)b*x.a)%p;
        a=(ll)a*x.a%p;
    }
}v;
vector<node> q[N];
int h[N],l[N],r[N],t;
void dfs(int x){
    if(!h[x]){
        q[r[x]=++t].push_back(v);
        return;
    }
    int m=0,y,z;
    for(int i=h[x];i;i=l[i]){
        dfs(i);
        if(i==h[x]){r[x]=r[i];continue;}
        y=q[r[x]].size(),z=q[r[i]].size();
        m=max(m,min(y,z));
        if(y<z) swap(y,z),swap(r[x],r[i]);
        for(int j=1;j<=z;++j)
            q[r[x]][y-j]+=q[r[i]][z-j];
    }
    m=q[r[x]].size()-m;
    for(int j=q[r[x]].size()-1;j>=m;--j){
        q[r[x]][j].a=(q[r[x]][j].a+q[r[x]][j].c)%p;
        q[r[x]][j].c=0;
    }
    q[r[x]].push_back(v);
}
ll P(int y){ll x=2,z=1;for(;y;y>>=1,x=x*x%p) if(y&1) z=z*x%p;return z;}
int main(){
    int n,x,s=0;
    scanf("%d",&n);
    v=(node){p+1>>1,p+1>>1,0};
    for(int i=1;i<=n;++i){
        scanf("%d",&x);
        l[i]=h[x],h[x]=i;
    }
    dfs(0);
    for(int i=q[r[0]].size()-1;i>=0;--i) s=(s+q[r[0]][i].b)%p;
    printf("%d\n",(ll)s*P(n+1)%p);
    return 0;
}