use proconio::input;

// modint {{{
#[allow(dead_code)]
mod modint {
    pub trait ModValue:
        std::ops::Add<Output = Self>
        + std::ops::Sub<Output = Self>
        + std::ops::Mul<Output = Self>
        + std::ops::Div<Output = Self>
        + std::ops::Rem<Output = Self>
        + std::ops::AddAssign
        + std::ops::SubAssign
        + std::ops::MulAssign
        + std::ops::DivAssign
        + std::ops::RemAssign
        + std::ops::Neg
        + std::cmp::Ord
        + Clone
        + Copy
        + Sized
        + std::fmt::Debug
        + std::fmt::Display
    {
    }

    impl<
            Value: std::ops::Add<Output = Self>
                + std::ops::Sub<Output = Self>
                + std::ops::Mul<Output = Self>
                + std::ops::Div<Output = Self>
                + std::ops::Rem<Output = Self>
                + std::ops::AddAssign
                + std::ops::SubAssign
                + std::ops::MulAssign
                + std::ops::DivAssign
                + std::ops::RemAssign
                + std::ops::Neg
                + std::cmp::Ord
                + Clone
                + Copy
                + Sized
                + std::fmt::Debug
                + std::fmt::Display,
        > ModValue for Value
    {
    }

    pub trait Minfo: Clone + Copy + std::fmt::Debug {
        type Value: ModValue;

        fn modulus() -> Self::Value;

        fn zero() -> Self::Value;

        fn one() -> Self::Value;
    }

    #[derive(Clone, Copy, Debug)]
    pub struct Mint<Mod: Minfo>(Mod::Value);

    impl<Mod: Minfo> Mint<Mod> {
        #[inline]
        fn is_within_the_range(x: Mod::Value) -> bool {
            Mod::zero() <= x && x < Mod::modulus()
        }

        #[inline]
        fn normalize(x: Mod::Value) -> Mod::Value {
            if Self::is_within_the_range(x) {
                x
            } else {
                let x = x % Mod::modulus();
                if x < Mod::zero() {
                    x + Mod::modulus()
                } else {
                    x
                }
            }
        }

        pub fn zero() -> Self {
            Self(Mod::zero())
        }

        pub fn one() -> Self {
            Self(Mod::one())
        }

        pub fn modulus() -> Mod::Value {
            Mod::modulus()
        }

        pub fn is_zero(&self) -> bool {
            self.value() == Mod::zero()
        }

        pub fn from_value(x: Mod::Value) -> Self {
            Self(Self::normalize(x))
        }

        pub fn from_frac(num: Mod::Value, den: Mod::Value) -> Self {
            Self::from_value(num) / Self::from_value(den)
        }

        pub fn value(&self) -> Mod::Value {
            self.0
        }

        pub fn inv(&self) -> Self {
            if self.is_zero() {
                panic!("inv(0)");
            }
            Self(Self::raw_inv(self.value()))
        }

        pub fn pow(&self, b: u64) -> Self {
            Self(Self::raw_pow(self.value(), b))
        }

        fn raw_add(x: Mod::Value, y: Mod::Value) -> Mod::Value {
            let z = x + y;
            if Mod::modulus() <= z {
                z - Mod::modulus()
            } else {
                z
            }
        }

        fn raw_sub(x: Mod::Value, y: Mod::Value) -> Mod::Value {
            let z = x - y;
            if z < Mod::zero() {
                z + Mod::modulus()
            } else {
                z
            }
        }

        fn raw_neg(x: Mod::Value) -> Mod::Value {
            if x == Mod::zero() {
                Mod::zero()
            } else {
                Mod::modulus() - x
            }
        }

        fn raw_mul(x: Mod::Value, y: Mod::Value) -> Mod::Value {
            x * y % Mod::modulus()
        }

        fn raw_div(x: Mod::Value, y: Mod::Value) -> Mod::Value {
            Self::raw_mul(x, Self::raw_inv(y))
        }

        #[allow(clippy::many_single_char_names)]
        fn raw_inv(mut x: Mod::Value) -> Mod::Value {
            let mut y = Mod::modulus();
            let mut u = Mod::one();
            let mut v = Mod::zero();
            while x != Mod::zero() {
                let q = y / x;
                y -= x * q;
                v -= u * q;
                std::mem::swap(&mut x, &mut y);
                std::mem::swap(&mut u, &mut v);
            }
            Self::normalize(v)
        }

        fn raw_pow(mut a: Mod::Value, mut b: u64) -> Mod::Value {
            let mut x = Mod::one();
            while 0 < b {
                if b % 2 == 1 {
                    x = Self::raw_mul(x, a);
                }
                a = Self::raw_mul(a, a);
                b /= 2
            }
            x
        }
    }

    macro_rules! impl_biop {
        (
            $(
                $biop_trait: ident::$biop_fn: ident,
                $biop_assign_trait: ident::$biop_assign_fn: ident,
                $biop_impl: ident;
            )*
        ) => {
            $(
                impl<Mod: Minfo> std::ops::$biop_trait for Mint<Mod> {
                    type Output = Self;
                    fn $biop_fn(self, rhs: Self) -> Self::Output {
                        Mint(Self::$biop_impl(self.value(), rhs.value()))
                    }
                }
                impl<Mod: Minfo> std::ops::$biop_assign_trait for Mint<Mod>
                where
                    Self: Copy
                {
                    fn $biop_assign_fn(&mut self, rhs: Self) {
                        use std::ops::$biop_trait;
                        *self = Mint::$biop_fn(*self, rhs)
                    }
                }
            )*
        };
    }

    impl_biop! {
        Add::add, AddAssign::add_assign, raw_add;
        Sub::sub, SubAssign::sub_assign, raw_sub;
        Mul::mul, MulAssign::mul_assign, raw_mul;
        Div::div, DivAssign::div_assign, raw_div;
    }

    impl<Mod: Minfo> std::ops::Neg for Mint<Mod> {
        type Output = Self;
        fn neg(self) -> Self::Output {
            Self(Self::raw_neg(self.value()))
        }
    }

    impl<Mod: Minfo> std::cmp::PartialEq for Mint<Mod> {
        fn eq(&self, other: &Self) -> bool {
            self.value().eq(&other.value())
        }
    }
    impl<Mod: Minfo> std::cmp::Eq for Mint<Mod> {}

    impl<Mod: Minfo> std::fmt::Display for Mint<Mod> {
        fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
            write!(f, "{}", self.value())
        }
    }

    impl<Mod: Minfo> std::iter::Sum<Self> for Mint<Mod> {
        fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
            iter.fold(Mint::zero(), |acc, x| acc + x)
        }
    }

    impl<Mod: Minfo> std::iter::Product<Self> for Mint<Mod> {
        fn product<I: Iterator<Item = Self>>(iter: I) -> Self {
            iter.fold(Mint::one(), |acc, x| acc * x)
        }
    }

    pub struct Factorial<Mod: Minfo> {
        normal: Vec<Mint<Mod>>,
        inverted: Vec<Mint<Mod>>,
    }

    impl<Mod: Minfo> Factorial<Mod>
    where
        Mod::Value: std::convert::TryFrom<usize, Error = std::num::TryFromIntError>,
    {
        pub fn new(n: usize) -> Self {
            if n == 0 {
                Self {
                    normal: Vec::new(),
                    inverted: Vec::new(),
                }
            } else {
                use std::convert::TryFrom;
                let mut normal = vec![Mint::one(); n];
                for i in 1..n {
                    normal[i] = normal[i - 1] * Mint::from_value(Mod::Value::try_from(i).unwrap());
                }
                let mut inverted = vec![Mint::one(); n];
                inverted[n - 1] = normal[n - 1].inv();
                for i in (1..n).rev() {
                    inverted[i - 1] =
                        inverted[i] * Mint::from_value(Mod::Value::try_from(i).unwrap());
                }
                Self { normal, inverted }
            }
        }
    }

    impl<Mod: Minfo> Factorial<Mod> {
        pub fn is_empty(&self) -> bool {
            self.normal.is_empty()
        }

        pub fn len(&self) -> usize {
            self.normal.len()
        }

        pub fn inverted(&self, i: usize) -> Mint<Mod> {
            self.inverted[i]
        }

        pub fn falling(&self, x: usize, p: usize) -> Mint<Mod> {
            if p == 0 {
                Mint::one()
            } else if x < p {
                Mint::zero()
            } else {
                assert!(x < self.len());
                self[x] * self.inverted(x - p)
            }
        }

        pub fn binom(&self, n: usize, k: usize) -> Mint<Mod> {
            if n < k {
                Mint::zero()
            } else {
                assert!(n < self.len());
                self[n] * self.inverted(k) * self.inverted(n - k)
            }
        }
    }

    impl<Mod: Minfo, I: std::slice::SliceIndex<[Mint<Mod>]>> std::ops::Index<I> for Factorial<Mod> {
        type Output = I::Output;

        #[inline]
        fn index(&self, index: I) -> &Self::Output {
            &self.normal[index]
        }
    }

    pub struct DoubleFactorial<Mod: Minfo> {
        normal: Vec<Mint<Mod>>,
        inverted: Vec<Mint<Mod>>,
    }

    impl<Mod: Minfo> DoubleFactorial<Mod>
    where
        Mod::Value: std::convert::TryFrom<usize, Error = std::num::TryFromIntError>,
    {
        pub fn new(n: usize) -> Self {
            use std::convert::TryFrom;
            let mut normal = vec![Mint::one(); n];
            for i in 2..n {
                normal[i] = normal[i - 2] * Mint::from_value(Mod::Value::try_from(i).unwrap());
            }
            let mut inverted = vec![Mint::one(); n];
            inverted[n - 1] = normal[n - 1].inv();
            inverted[n - 2] = normal[n - 2].inv();
            for i in (2..n).rev() {
                inverted[i - 2] = inverted[i] * Mint::from_value(Mod::Value::try_from(i).unwrap());
            }
            Self { normal, inverted }
        }
    }

    impl<Mod: Minfo> DoubleFactorial<Mod> {
        pub fn inverted(&self, i: usize) -> Mint<Mod> {
            self.inverted[i]
        }
    }

    impl<Mod: Minfo, I: std::slice::SliceIndex<[Mint<Mod>]>> std::ops::Index<I>
        for DoubleFactorial<Mod>
    {
        type Output = I::Output;

        #[inline]
        fn index(&self, index: I) -> &Self::Output {
            &self.normal[index]
        }
    }

    #[macro_export]
    macro_rules! define_mint {
        (
            struct $modint: ident ($modinfo: ident ($modvalue: expr; $modty: ty));
        ) => {
            #[derive(Clone, Copy, Debug)]
            pub struct $modinfo {}
            impl Minfo for $modinfo {
                type Value = $modty;
                fn modulus() -> Self::Value {
                    $modvalue
                }
                fn zero() -> Self::Value {
                    0
                }
                fn one() -> Self::Value {
                    1
                }
            }
            type $modint = Mint<$modinfo>;
        };
    }

    #[derive(Clone, Copy, Debug)]
    pub struct Minfo1000000007 {}

    impl Minfo for Minfo1000000007 {
        type Value = i64;
        fn modulus() -> Self::Value {
            1_000_000_007
        }
        fn zero() -> Self::Value {
            0
        }
        fn one() -> Self::Value {
            1
        }
    }

    pub type Mint1000000007 = Mint<Minfo1000000007>;

    #[derive(Clone, Copy, Debug)]
    pub struct Minfo998244353 {}

    impl Minfo for Minfo998244353 {
        type Value = i64;
        fn modulus() -> Self::Value {
            998_244_353
        }
        fn zero() -> Self::Value {
            0
        }
        fn one() -> Self::Value {
            1
        }
    }

    pub type Mint998244353 = Mint<Minfo998244353>;
}
// }}}

type Mint = modint::Mint1000000007;

fn cal(a: &[Mint]) -> Mint {
    let n = a.len();
    let mut left = vec![Mint::one(); n + 1];
    let mut right = vec![Mint::one(); n + 1];
    for i in 0..n {
        left[i + 1] = left[i] * (Mint::one() - a[i]);
    }
    for i in 0..n {
        right[i + 1] = right[i] * (Mint::one() - a[n - 1 - i]);
    }
    a.iter()
        .enumerate()
        .map(|(i, &x)| left[i] * x * right[n - i - 1])
        .sum::<Mint>()
}

fn merge(first: &mut Vec<Mint>, f: &[Vec<Mint>]) {
    if !f.is_empty() {
        let second = &f[0];
        for i in 0..second.len() {
            let values: Vec<Mint> = f
                .iter()
                .take_while(|v| (second.len() - v.len()) <= i)
                .map(|v| v[i - (second.len() - v.len())])
                .chain(std::iter::once(first[i + first.len() - second.len()]))
                .collect();
            let j = i + first.len() - second.len();
            first[j] = cal(&values);
        }
    }
    first.push(Mint::from_value(2).inv());
}

fn main() {
    input!(n: usize, p: [usize; n]);
    let n = n + 1;
    let mut g = vec![vec![]; n + 1];
    for (i, x) in p.into_iter().enumerate().map(|(i, x)| (i + 1, x)) {
        g[x].push(i);
    }

    fn dfs(x: usize, g: &[Vec<usize>]) -> Vec<Mint> {
        let mut f: Vec<Vec<Mint>> = g[x].iter().map(|&y| dfs(y, &g)).collect();
        if f.is_empty() {
            vec![Mint::from_value(2).inv()]
        } else {
            f.sort_by_key(|v| std::cmp::Reverse(v.len()));
            let (first, others) = f.split_first_mut().unwrap();
            merge(first, others);
            let mut ans = Vec::new();
            std::mem::swap(&mut ans, &mut f[0]);
            ans
        }
    }
    let ans = dfs(0, &g);
    println!(
        "{}",
        ans.into_iter().sum::<Mint>() * Mint::from_value(2).pow(n as u64)
    );
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_cal_one() {
        assert_eq!(cal(&[Mint::from_value(42)]), Mint::from_value(42));
    }

    #[test]
    fn test_cal_two() {
        for &(a, b) in [(3, 4), (345, 23525)].iter() {
            let a = Mint::from_value(a);
            let b = Mint::from_value(b);
            assert_eq!(cal(&[a, b]), a * (Mint::one() - b) + (Mint::one() - a) * b);
        }
    }

    #[test]
    fn test_cal_three() {
        for &(a, b, c) in [(3, 4, 3), (13345, 23525, 2224)].iter() {
            let a = Mint::from_value(a);
            let b = Mint::from_value(b);
            let c = Mint::from_value(c);
            assert_eq!(
                cal(&[a, b, c]),
                a * (Mint::one() - b) * (Mint::one() - c)
                    + (Mint::one() - a) * b * (Mint::one() - c)
                    + (Mint::one() - a) * (Mint::one() - b) * c
            );
        }
    }

    #[test]
    fn test_merge_one() {
        for a in [vec![1, 3, 4], vec![534, 4]]
            .iter()
            .map(|v| v.iter().copied().map(Mint::from_value).collect::<Vec<_>>())
        {
            let mut first = a.clone();
            merge(&mut first, &Vec::new());
            assert_eq!(
                first,
                a.iter()
                    .copied()
                    .chain(std::iter::once(Mint::from_value(2).inv()))
                    .collect::<Vec<_>>()
            );
        }
    }

    #[test]
    fn test_merge_two() {
        let a = vec![Mint::from_value(2).inv()];
        let mut first = a.clone();
        merge(&mut first, &[vec![Mint::from_value(2).inv()]]);
        assert_eq!(
            first,
            vec![Mint::from_value(2).inv(), Mint::from_value(2).inv()]
        );
    }

    #[test]
    fn ters_merge_sample_3_vertex_4() {
        let half = Mint::from_frac(1, 2);
        let a = vec![half, half, half, half];
        let mut first = a.clone();
        merge(
            &mut first,
            &[vec![half, half, half], vec![half, half, half]],
        );
        assert_eq!(
            first,
            vec![
                half,
                Mint::from_frac(3, 8),
                Mint::from_frac(3, 8),
                Mint::from_frac(3, 8),
                half
            ]
        );
    }
}


Measured memory_bytes: 88164