#include<bits/stdc++.h>
using namespace std;
const int MAXN=5e5,MOD=998244353,g=3;
#define ll long long
int rev[MAXN<<1];
ll w[(1<<21)+5],ninv[MAXN],fac[MAXN],ifac[MAXN];
ll qpow(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1)ans*=a,ans%=MOD;
		a*=a,a%=MOD;b>>=1;
	}
	return ans;
}
void init()
{
	ll stp=qpow(g,(MOD-1)/(1<<21));
	w[0]=1;
	for(int i=1;i<(1<<21);i++)
		w[i]=w[i-1]*stp%MOD;
	ninv[1]=1;
	for(int i=2;i<=MAXN-5;i++)
		ninv[i]=(MOD-MOD/i)*ninv[MOD%i]%MOD;
	fac[0]=ifac[0]=1;
	for(int i=1;i<=MAXN-5;i++)
		fac[i]=fac[i-1]*i%MOD,
		ifac[i]=ifac[i-1]*ninv[i]%MOD;
}
struct node
{
	int len;
	vector<ll> a;
	void rsiz(int x)
	{
		len=x;
		a.resize(x);
	}
	void srk()
	{
		for(int i=len-1;i>=0;i--)
			if(a[i])
				{len=i+1;break;}
		a.resize(len);
	}
};
node makenode(int c)
{
	node x;x.rsiz(1);
	x.a[0]=c;
	return x;
}
node operator+(node x,node y)
{
	if(x.len<y.len)swap(x,y);
	for(int i=0;i<y.len;++i)
	{
		x.a[i]+=y.a[i];
		if(x.a[i]>MOD)x.a[i]-=MOD;
	}
	return x;
}
node operator-(node x,node y)
{
	int flag=0;
	if(x.len<y.len)swap(x,y),flag=1;
	for(int i=0;i<y.len;i++)
	{
		if(flag)x.a[i]=y.a[i]-x.a[i];
		if(x.a[i]<0)x.a[i]+=MOD;
	}
	for(int i=y.len;i<x.len;i++)
		x.a[i]=MOD-x.a[i];
	return x;
}
node operator*(node x,ll c)
{
	for(int i=0;i<x.len;i++)
		x.a[i]=x.a[i]*c%MOD;
	return x;
}
node operator*(ll c,node x)
{
	for(int i=0;i<x.len;i++)
		x.a[i]=x.a[i]*c%MOD;
	return x;
}
node operator>>(node x,int c)
{
	for(int i=0;i<=x.len-c-1;i++)
		x.a[i]=x.a[i+c];
	return x;
}
void pre(int rev[],int len)
{
	rev[0]=0;
	rev[1]=(len>>1);
	for(int i=2;i<len;i++)
	{
		rev[i]=rev[i>>1]>>1;
		if(i&1)rev[i]+=(len>>1);
	}
}
void ntt(node &x,int r)
{
	int len=x.len;
	pre(rev,len);
	for(int i=0;i<len;i++)
		if(i<rev[i])swap(x.a[i],x.a[rev[i]]);
	for(int k=2;k<=len;k<<=1)
	{
		int d=k>>1,step=(1<<21)/k;
		for(int i=0;i<len;i+=k)
		{
			int wn=0;
			for(int j=i;j<i+d;j++)
			{
				ll u=x.a[j],v=x.a[j+d];
				x.a[j]=(u+w[wn]*v)%MOD;
				x.a[j+d]=(u-w[wn]*v)%MOD;
				wn+=step;
			}
		}
	}
	if(r==-1)
	{
		reverse(x.a.begin()+1,x.a.end());
		ll inv=qpow(len,MOD-2);
		for(int i=0;i<len;i++)
			x.a[i]=x.a[i]*inv%MOD;
	}
}
int getlen(int x,int y)
{
	x+=y-1;
	int res=1;
	while(res<x)res<<=1;
	return res;
}
node operator*(node x,node y)
{
	int len=getlen(x.len,y.len);
	x.rsiz(len),y.rsiz(len);
	ntt(x,1),ntt(y,1);
	for(int i=0;i<len;i++)
		x.a[i]=x.a[i]*y.a[i]%MOD;
	ntt(x,-1);
	x.srk();
	return x;
}
node inv(node x)
{
	if(x.len==1)
	{
		x.a[0]=qpow(x.a[0],MOD-2);
		return x;
	}
	node y=x;y.rsiz((x.len+1)/2);
	node f0=inv(y);
	int xlen=x.len,len=getlen(f0.len*2,x.len);
	f0.rsiz(len),x.rsiz(len);
	ntt(f0,1),ntt(x,1);
	for(int i=0;i<len;i++)
		x.a[i]=(2*f0.a[i]-f0.a[i]*f0.a[i]%MOD*x.a[i])%MOD;
	ntt(x,-1);
	x.rsiz(xlen);
	return x;
}
node sqr(node x)
{
	if(x.len==1)
	{
		x.a[0]=1;
		return x;
	}
	node y=x;y.rsiz((x.len+1)/2);
	node f0=sqr(y);
	int xlen=x.len,len=getlen(xlen,f0.len*2);
	f0.rsiz(len);
	node f1=inv(f0);
	len=getlen(len,f1.len);
	x.rsiz(len);f1.rsiz(len);f0.rsiz(len);
	ntt(f0,1),ntt(f1,1),ntt(x,1);
	ll hlf=qpow(2,MOD-2);
	for(int i=0;i<len;i++)
		x.a[i]=(f0.a[i]*f0.a[i]%MOD+x.a[i])*hlf%MOD,
		x.a[i]=x.a[i]*f1.a[i]%MOD;
	ntt(x,-1);
	x.rsiz(xlen);
	return x;
}
node drv(node x)
{
	int len=x.len;
	for(int i=0;i<len-1;i++)
		x.a[i]=x.a[i+1]*(i+1)%MOD;
	x.rsiz(len-1);
	return x;
}
node idrv(node x)
{
	int len=x.len;
	x.rsiz(len+1);
	for(int i=len-1;i>=1;i--)
		x.a[i]=x.a[i-1]*qpow(i,MOD-2)%MOD;
	x.a[0]=0;
	return x;
}
node lne(node x)
{
	int len=x.len;
	x=drv(x)*inv(x);
	x.rsiz(len);x=idrv(x);
	x.rsiz(len);
	return x;
}
node exp(node x)
{
	if(x.len==1)
	{
		x.a[0]=1;
		return x;
	}
	int len=x.len;
	node y=x;y.rsiz((len+1)/2);
	node f0=exp(y);f0.rsiz(len);
	node f1=lne(f0);
	x=f0*(makenode(1)-f1+x);x.rsiz(len);
	return x;
}
int n,fa[MAXN],dep[MAXN];
node f[MAXN];
vector<int> ve[MAXN];
void dfs(int o)
{
	dep[o]=dep[fa[o]]+1;
	for(auto u:ve[o])
		dfs(u);
}
void calc(int o)
{
	f[o].rsiz(dep[o]+1);
	for(int i=0;i<=dep[o];i++)
		f[o].a[i]=ifac[i];
	for(auto u:ve[o])
	{
		calc(u);
		f[o]=f[o]*f[u];
	}
	f[o]=f[o]>>1;
	f[o].rsiz(dep[o]);
}
int main()
{
	init();
	cin>>n;
	for(int i=2;i<=n;i++)
		scanf("%d",fa+i),ve[fa[i]].push_back(i);
	dfs(1);
	calc(1);
	cout<<(f[1].a[0]*fac[n]%MOD+MOD)%MOD;
	return 0;
}