Snuke has a rooted tree with N+1 vertices.
The vertices are numbered 0 through N, and Vertex 0 is the root of the tree.
The parent of Vertex i (1 \leq i \leq N) is Vertex p_i.
Besides this tree, Snuke also has an box which is initially empty and many marbles, and playing with them.
The play begins with placing one marble on some of the vertices, then proceeds as follows:
 - If there is a marble on Vertex 0, move the marble into the box.
 - Move each marble from the vertex to its parent (all at once).
 - For each vertex occupied by two or more marbles, remove all the marbles from the vertex.
 - If there exists a vertex with some marbles, go to Step 1. Otherwise, end the play.
There are 2^{N+1} ways to place marbles on some of the vertices.
For each of them, find the number of marbles that will be in the box at the end of the play, and compute the sum of all those numbers modulo 1,000,000,007.

-----Constraints-----
 - 1 \leq N < 2 \times 10^{5}
 - 0 \leq p_i < i

-----Partial Scores-----
 - In the test set worth 400 points, N < 2{,}000.

-----Input-----
Input is given from Standard Input in the following format:
N
p_1 p_2 ... p_{N}

-----Output-----
Print the answer.

-----Sample Input-----
2
0 0

-----Sample Output-----
8

When we place a marble on both Vertex 1 and 2, there will be multiple marbles on Vertex 0 by step 2. In such a case, these marbles will be removed instead of being moved to the box.

Difficulty: competition

Starter code:

Solutions:
# coding: utf-8
# Your code here!
import sys
read = sys.stdin.read
readline = sys.stdin.readline

n, = list(map(int, readline().split()))
p = [-1] + [*list(map(int, readline().split()))]

MOD = 10**9+7
child = [[] for i in range(n+1)]
tot = [None for i in range(n+1)]
one = [None for i in range(n+1)]
dep = [0]*(n+1)
p2 = [1]*(n+1)
for i in range(n):
    p2[i+1] = p2[i]*2%MOD

for v in range(n,-1,-1):
    if dep[v]==0:
        tot[v] = []
        one[v] = []
    else:
        child[v].sort(key=lambda i: dep[i])
        one[v] = one[child[v][-1]]
        tot[v] = tot[child[v][-1]]
        #one_sum = [0]*(dep[v])
        #zero_sum = [0]*(dep[v])
        child[v].pop()
        if child[v]:
            zero = [p2[tot[v][j]]-one[v][j] for j in range(-len(one[child[v][-1]]),0)]
        for c in child[v]:
            for j in range(-len(one[c]),0):
                z = p2[tot[c][j]]-one[c][j]
                one[v][j] = (one[v][j]*z+zero[j]*one[c][j])%MOD
                zero[j] = zero[j]*z%MOD
                tot[v][j] += tot[c][j]

    tot[v].append(1)
    one[v].append(1)

    child[p[v]].append(v)
    dep[p[v]] = max(dep[p[v]],dep[v]+1)        

    #print(v,tot[v],one[v])
        
#print("tot",tot[0])
#print("one",one[0])

ans = 0
for i,j in zip(tot[0],one[0]):
    ans += pow(2,n+1-i,MOD)*j%MOD

print((ans%MOD))
#print(sum(tot[0]))



"""

https://atcoder.jp/contests/arc086/tasks/arc086_c

対消滅する可能性があるのは同じ深さの点に置かれたビー玉だけ
→1つも残らない or 1つだけ残るである

すなわち、ある深さに関して、そのうち1つ残るのがいくつあるかを数えればよい

1つだけ置く場合→絶対のこる
2つ置く場合→絶対消える
3つ置く場合→3つのLCAが等しくなければ残る

木dpみたいにする？
dp[0] = 一つも含んでいない場合の通り数
dp[1] = 1つ含んでいる場合の通り数

dp[0] = すべての場合 - dp[1]
dp[1] = 1つだけ1を持っている場合

で、子から親に伝播させていく…？
マージテクで計算量削減か？
3sだからそうっぽいな…
3つ以上のマージ書くのだるすぎん？

またはマージの順番をどっかにメモっておく

maxd-d = indexでやる
最長のやつにマージする

dp[maxd-d][1] = 1つだけ元深さdが残っている場合の数
c = 子の数

計算量は？
サイズは高々1しか増えないので可能っぽい
2倍処理がまずい

もっと簡潔に？
単一のdだけで考えよう
2倍処理なんてしない
最後に各dに関してかければいい
→するとマージの際に2番目の大きさだけでなんとかなる

必要のないマージをしない
両方に関係ないdは操作しない
むしろばらつきは深い部分にだけ存在するか
浅い部分は共通。よって-xで管理すればいいか

"""

import sys
mod = 10**9 + 7
sys.setrecursionlimit(200000)

from collections import deque
def NC_Dij(lis,start):

    ret = [float("inf")] * len(lis)
    ret[start] = 0
    
    q = deque([start])
    plis = [i for i in range(len(lis))]

    while len(q) > 0:
        now = q.popleft()

        for nex in lis[now]:

            if ret[nex] > ret[now] + 1:
                ret[nex] = ret[now] + 1
                plis[nex] = now
                q.append(nex)

    return ret,plis

def inverse(a): #aのmodを法にした逆元を返す
    return pow(a,mod-2,mod)


def dfs(v):

    if len(lis[v]) == 0:
        ret = [ [1,1] ]
        return ret

    else:

        retlis = []
        for nex in lis[v]:
            nret = dfs(nex)
            retlis.append( [len(nret),nret] )
        retlis.sort()

        #1つしかない場合マージしない
        if len(retlis) == 1:
            retlis[-1][1].append([1,1])
            return retlis[-1][1]

        #2つ以上の場合最大のやつにマージする
        for revd in range(retlis[-2][0]):

            zmul = 1
            amul = 1
            for i in range(len(retlis)-1,-1,-1):
                if revd < retlis[i][0]:
                    zmul *= retlis[i][1][-1-revd][0]
                    amul *= sum(retlis[i][1][-1-revd])
                    zmul %= mod
                    amul %= mod
                else:
                    break

            nsum = 0
            for i in range(len(retlis)-1,-1,-1):
                if revd < retlis[i][0]:
                    nsum += zmul * inverse(retlis[i][1][-1-revd][0]) * retlis[i][1][-1-revd][1]
                    nsum %= mod
                else:
                    break

            retlis[-1][1][-1-revd][1] = nsum
            retlis[-1][1][-1-revd][0] = (amul-nsum) % mod

        retlis[-1][1].append([1,1])
        return retlis[-1][1]

N = int(input())
p = list(map(int,input().split()))

lis = [ [] for i in range(N+1) ]

for i in range(N):

    #lis[i+1].append(p[i])
    lis[p[i]].append(i+1)

dlis,plis = NC_Dij(lis,0)
maxd = max(dlis)

dn = [0] * (maxd+1)
for i in dlis:
    dn[i] += 1

ans = dfs(0)
#print (dn,ans)
A = 0
for i in range(maxd+1):
    A += ans[-1-i][1] * pow(2,N+1-dn[i],mod)
    A %= mod
print (A)


import sys
readline = sys.stdin.readline

 
def parorder(Edge, p):
    N = len(Edge)
    par = [0]*N
    par[p] = -1
    stack = [p]
    order = []
    visited = set([p])
    ast = stack.append
    apo = order.append
    while stack:
        vn = stack.pop()
        apo(vn)
        for vf in Edge[vn]:
            if vf in visited:
                continue
            visited.add(vf)
            par[vf] = vn
            ast(vf)
    return par, order
 
def getcld(p):
    res = [[] for _ in range(len(p))]
    for i, v in enumerate(p[1:], 1):
        res[v].append(i)
    return res
 
 
 
MOD = 10**9+7
limit = 1341398
p2 = [1]*limit
for i in range(1, limit):
    p2[i] = 2*p2[i-1]%MOD

N = int(readline())
P = [None] + list(map(int, readline().split()))
Edge = [[] for _ in range(N+1)]
for i in range(1, N+1):
    Edge[i].append(P[i])
    Edge[P[i]].append(i)
_, L = parorder(Edge, 0)
C = getcld(P)
dist = [0]*(N+1)
for l in L[1:]:
    p = P[l]
    dist[l] = 1 + dist[p]
md = dist[:]
for l in L[:0:-1]:
    p = P[l]
    md[p] = max(md[p], md[l])

Ci = [None]*(N+1)
for i in range(N+1):
    if C[i]:
        res = -1
        cc = None
        for ci in C[i]:
            if md[ci] > res:
                res = md[ci]
                cc = ci
    Ci[i] = cc
ans = 0

res = [[0]]
val = [1]*(N+1) 
size = [1]*(N+1) 
used = [0]*(N+1)
name = [None]*(N+1)
mark = [None]*(N+1)
while res:
    newres = []
    for a in res:
        ts = 1
        ks = 0
        ss = 0
        dic = []
        ctr = []

        for i in a:
            ts = ts*(p2[size[i]] - val[i])%MOD
            ks = (ks + val[i]*pow(p2[size[i]] - val[i], MOD-2, MOD))%MOD
            ss += size[i]

            if not used[i]:
                used[i] = 1
                if C[i]:
                    dic.append(C[i])
                    if name[i] is None:
                        name[i] = Ci[i]
            if name[i] is not None:
                ctr.append(name[i])
        newres.extend(dic)
        if len(ctr) > 1:
            newres.append(ctr)
        for a0 in a:
            val[a0] = ts*ks%MOD
            size[a0] = ss
    ans = (ans + val[a0]*p2[N+1 - size[a0]])%MOD
    res = [nr[:] for nr in newres]
    #print(res)
print(ans)

"""
def merge(sz, vl):
    res = 0
    ts = 1
    ks = 0
    for s, v in zip(sz, vl):
        ts = ts*(p2[s] - v)%MOD
        ks = (ks + v*pow(p2[s] - v, MOD-2, MOD))%MOD
    return ts*ks%MOD



assert N < 2000

st = [0]
aa = 0
while st:
    size = [0]*(N+1)
    val = [0]*(N+1)
    cs = [[] for _ in range(N+1)]
    cv = [[] for _ in range(N+1)]
    
    nr = []
    for s in st:
        cs[s] = [1]
        cv[s] = [1]
        nr.extend(C[s])
    for l in L[::-1]:
        size[l] = sum(cs[l])
        val[l] = merge(cs[l], cv[l])
        if l:
            p = P[l]
            cs[p].append(size[l])
            cv[p].append(val[l])
    aa = (aa + val[0]*p2[N+1 - size[0]])%MOD
    st = nr[:]

print(aa)
"""
    


from collections import deque


def get_pow():
    cache = {}

    def func(x):
        if x not in cache:
            cache[x] = pow(2, x, mod)
        return cache[x]

    return func


mod = 1000000007
n = int(input())
parents = list(map(int, input().split()))
children = [set() for _ in range(n + 1)]
for c, p in enumerate(parents):
    children[p].add(c + 1)

levels = [{0}]
while True:
    level = set()
    for p in levels[-1]:
        level.update(children[p])
    if not level:
        break
    levels.append(level)
levels.reverse()

level_node_count = []
balls = [None] * (n + 1)

for i, level in enumerate(levels):
    level_node_count.append(len(level))
    for node in level:
        cn = children[node]
        if cn:
            if len(cn) == 1:
                bs = balls[cn.pop()]
                bs.appendleft([1, 1, 0])
                balls[node] = bs
                continue
            balls_from_children = [balls[c] for c in children[node]]
            balls_from_children.sort(key=len)
            bs1 = balls_from_children[0]
            for bs2 in balls_from_children[1:]:
                for (b10, b11, b12), b2 in zip(bs1, bs2):
                    b2[2] = ((b11 + b12) * b2[1] + b12 * b2[0]) % mod
                    b2[1] = (b10 * b2[1] + b11 * b2[0]) % mod
                    b2[0] = b2[0] * b10 % mod
                bs1 = bs2
            for b in bs1:
                b[0] = (b[0] + b[2]) % mod
                b[2] = 0
            bs1.appendleft([1, 1, 0])
            balls[node] = bs1
        else:
            balls[node] = deque([[1, 1, 0]])

level_node_count.reverse()

pow2 = get_pow()
print((sum(b[1] * pow2(n - l + 1) % mod for l, b in zip(level_node_count, balls[0])) % mod))


from collections import deque


def get_pow():
    cache = {}

    def func(x):
        if x not in cache:
            cache[x] = pow(2, x, mod)
        return cache[x]

    return func


mod = 1000000007
n = int(input())
parents = list(map(int, input().split()))
children = [set() for _ in range(n + 1)]
for c, p in enumerate(parents):
    children[p].add(c + 1)

levels = [{0}]
while True:
    level = set()
    for p in levels[-1]:
        level.update(children[p])
    if not level:
        break
    levels.append(level)
levels.reverse()

level_node_count = []
balls = [None] * (n + 1)

for i, level in enumerate(levels):
    level_node_count.append(len(level))
    for node in level:
        cn = children[node]
        if cn:
            if len(cn) == 1:
                bs = balls[cn.pop()]
                bs.appendleft([1, 1, 0])
                balls[node] = bs
                continue
            balls_from_children = [balls[c] for c in children[node]]
            balls_from_children.sort(key=len)
            bs1 = balls_from_children[0]
            for bs2 in balls_from_children[1:]:
                for (b10, b11, b12), b2 in zip(bs1, bs2):
                    b2[2] = ((b11 + b12) * b2[1] + b12 * b2[0]) % mod
                    b2[1] = (b10 * b2[1] + b11 * b2[0]) % mod
                    b2[0] = b2[0] * b10 % mod
                bs1 = bs2
            lim = len(balls_from_children[-2])
            for i, b in enumerate(bs1):
                if i >= lim:
                    break
                b[0] = (b[0] + b[2]) % mod
                b[2] = 0
            bs1.appendleft([1, 1, 0])
            balls[node] = bs1
        else:
            balls[node] = deque([[1, 1, 0]])

level_node_count.reverse()

pow2 = get_pow()
print((sum(b[1] * pow2(n - l + 1) % mod for l, b in zip(level_node_count, balls[0])) % mod))


# seishin.py
from collections import deque
N = int(input())
*P, = list(map(int, input().split()))
MOD = 10**9 + 7

G = [[] for i in range(N+1)]
U = [0]*(N+1)
C = [0]*(N+1)
for i, p in enumerate(P):
    G[p].append(i+1)
    U[i+1] = u = U[p]+1
    C[u] += 1
Q = [None]*(N+1)
PP = {}
def pp(k):
    if k not in PP:
        PP[k] = p = pow(2, k, MOD)
        return p
    return PP[k]
L = [0]*(N+1)

ept = []
sz = L.__getitem__
for i in range(N, -1, -1):
    g = G[i]
    if not g:
        continue

    # 子ノードのdequeを集める
    g.sort(key=sz, reverse=1)
    k = len(g)
    e = [pp(k) - k, k, 0]

    g0 = g[0]
    L[i] = L[g0] + 1
    if L[g0] == 0:
        Q[i] = deque([e])
        continue
    Q[i] = R = Q[g0]
    if k > 1:
        # a0 <- a2
        for s, r in zip(Q[g[1]] or ept, R):
            r[0] += r[2]; r[2] = 0

        for j in g[1:]:
            S = Q[j]
            if not S:
                break

            # dequeの小さい方から大きい方へマージする処理
            for (a0, a1, a2), r in zip(S, R):
                b0, b1, b2 = r; a0 += a2
                r[0] = a0*b0 % MOD
                r[1] = (a0*b1 + a1*b0) % MOD
                r[2] = ((a0+a1)*b2 + a1*b1) % MOD
    R.appendleft(e)
print(((pp(N) + sum(pp(N+1-c) * a1 % MOD for (a0, a1, a2), c in zip(Q[0], C[1:]))) % MOD))


import sys
input = sys.stdin.readline
sys.setrecursionlimit(10**7)
from collections import deque

MOD = 10 ** 9 + 7

N = int(input())
graph = [[] for _ in range(N+1)]
for i,x in enumerate(input().rstrip().split(),1):
    x = int(x)
    graph[i].append(x)
    graph[x].append(i)

# 各深さから来ている (0個、1個、2個以上) の分布を 確率 mod MODで持つ。

half = (MOD + 1) // 2

def merge(dp,dp1):
    L = len(dp1)
    for i in range(L):
        # 0個,1個,2個以上
        a,b,c = dp[i]
        d,e,f = dp1[i]
        a,b,c = a*d, a*e + b*d, a*f + b*e + b*f + c*d + c*e + c*f
        a %= MOD
        b %= MOD
        c %= MOD
        dp[i] = (a,b,c)
    return

def dfs(v,parent = None):
    dp = None
    L = 0
    for u in graph[v]:
        if u == parent:
            continue
        dp1 = dfs(u,v)
        if dp is None:
            dp = dp1
        else:
            if len(dp) < len(dp1):
                dp,dp1 = dp1,dp
            # 2個以上が入っているインデックス
            if L < len(dp1):
                L = len(dp1)
            merge(dp,dp1)
    if dp is None:
        dp = deque()
    else:
        # 2個以上あるときに、0個化する
        for i in range(L):
            a,b,c = dp[i]
            dp[i] = (a+c,b,0)
    dp.appendleft((half,half,0))
    return dp

dp = dfs(0)

answer = sum(b for a,b,c in dp)
answer *= pow(2,N+1,MOD)
answer %= MOD
print(answer)

import sys

sys.setrecursionlimit(10 ** 6)
input = sys.stdin.readline
int1 = lambda x: int(x) - 1
p2D = lambda x: print(*x, sep="\n")

def main():
    def dfs(u=0):
        def merge(dpu, dpv):
            vn = len(dpv)
            for d in range(-1, -1 - vn, -1):
                u0, u1, u2 = dpu[d]
                v0, v1, v2 = dpv[d]
                n0 = (u0 * v0) % md
                n1 = (u0 * v1 + u1 * v0) % md
                n2 = (u2 * (v0 + v1 + v2) + v2 * (u0 + u1) + u1 * v1) % md
                dpu[d] = (n0, n1, n2)

        # 葉の場合
        if len(to[u]) == 0:
            return [(inv2, inv2, 0)]
        # すべての子をマージ
        dpu = []
        mxlen=0
        for v in to[u]:
            dpv = dfs(v)
            #深さが2段以上あったらu2をu0に
            if not dpu:
                dpu = dpv
            else:
                if len(dpu) < len(dpv): dpu, dpv = dpv, dpu
                mxlen=max(mxlen,len(dpv))
                merge(dpu, dpv)
        for d in range(-1,-1-mxlen,-1):
            u0,u1,u2=dpu[d]
            dpu[d] = (u0 + u2, u1, 0)
        dpu.append((inv2, inv2, 0))
        return dpu

    md = 10 ** 9 + 7
    # 1/2のmod
    inv2 = pow(2, md - 2, md)
    n = int(input())
    to = [[] for _ in range(n+1)]
    pp = list(map(int, input().split()))
    for i, p in enumerate(pp, 1):
        to[p].append(i)
    # print(to)
    dp0 = dfs()
    # print(dp0)
    ans = sum(u1 for _, u1, _ in dp0)
    print((ans * pow(2, n + 1, md)) % md)

main()


n, = map(int, input().split())
p = [-1] + [*map(int, input().split())]

MOD = 10**9+7
dp = [[] for _ in range(n+1)]
dep = [0]*(n+1)
nxt = [0]*(n+1)

for v in range(n,0,-1):
    _,nxt[p[v]],dep[p[v]] = sorted([nxt[p[v]],dep[p[v]],dep[v]+1])

tot = [0]*(dep[0]+1)
for i in range(n+1): tot[dep[i]] += 1

def merge(p,v):
    if len(dp[p]) < len(dp[v]):
        dp[p],dp[v]=dp[v],dp[p]
    for i in range(-len(dp[v]),0):
        a,b,c = dp[p][i]
        d,e,f = dp[v][i]
        dp[p][i][:] = [a*d%MOD,(b*d+a*e)%MOD,c*f%MOD]

for v in range(n,-1,-1):
    dp[v].append([1,1,2])
    for i in range(-nxt[v]-1,0):
        dp[v][i][0] = dp[v][i][2] - dp[v][i][1]
    if v: merge(p[v],v)

ans = 0
for d in dp[0]:
    ans += pow(d[2],MOD-2,MOD)*d[1]%MOD
print(ans*pow(2,n+1,MOD)%MOD)

# seishin.py
from collections import deque
N = int(input())
*P, = map(int, input().split())
MOD = 10**9 + 7

G = [[] for i in range(N+1)]
U = [0]*(N+1)
C = [0]*(N+1)
for i, p in enumerate(P):
    G[p].append(i+1)
    U[i+1] = u = U[p]+1
    C[u] += 1
Q = [None]*(N+1)
PP = {}
def pp(k):
    if k not in PP:
        PP[k] = p = pow(2, k, MOD)
        return p
    return PP[k]
L = [0]*(N+1)

ept = []
sz = L.__getitem__
for i in range(N, -1, -1):
    g = G[i]
    if not g:
        continue

    # 子ノードのdequeを集める
    g.sort(key=sz, reverse=1)
    k = len(g)
    e = [pp(k) - k, k, 0]

    g0 = g[0]
    L[i] = L[g0] + 1
    if L[g0] == 0:
        Q[i] = deque([e])
        continue
    Q[i] = R = Q[g0]
    if k > 1:
        # a0 <- a2
        for s, r in zip(Q[g[1]] or ept, R):
            r[0] += r[2]; r[2] = 0

        for j in g[1:]:
            S = Q[j]
            if not S:
                break

            # dequeの小さい方から大きい方へマージする処理
            for (a0, a1, a2), r in zip(S, R):
                b0, b1, b2 = r; a0 += a2
                r[0] = a0*b0 % MOD
                r[1] = (a0*b1 + a1*b0) % MOD
                r[2] = ((a0+a1)*b2 + a1*b1) % MOD
    R.appendleft(e)
print((pp(N) + sum(pp(N+1-c) * a1 % MOD for (a0, a1, a2), c in zip(Q[0], C[1:]))) % MOD)