#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define MP make_pair
mt19937 rnd(time(0));
#define ll long long
#define LL __int128
#define MP make_pair
#define poly vector<long long>
const int MAXN=5e5+5;
const int MOD=998244353;
const ll inv2=(MOD+1)/2;
namespace polynomial{// yjl poly plank 5.0 ver
int bfly[MAXN];ll inver[MAXN];
int clogg(int x){return (int)ceil(log2(x));} 
ll ksm(ll a,int b){ll res=1;while(b){if(b&1)res=res*a%MOD;a=a*a%MOD,b>>=1;}return res;}
void butterfly(int l){
	static int las;
	if(las!=l){
		las=l; 
		for(int i=1;i<(1<<l);i++)
			bfly[i]=(bfly[i>>1]>>1)|((i&1)<<l-1);
	} 
}
void NTT(poly &f,int l,int typ){
	butterfly(l);f.resize(1<<l);
	for(int i=0;i<(1<<l);i++)
		if(bfly[i]<i) swap(f[i],f[bfly[i]]);
	for(int i=0;i<l;i++){
		ll step=ksm(3,MOD-1+(MOD-1>>i+1)*typ);
		for(int j=0;j<(1<<l);j+=(1<<i+1)){
			ll cur=1,l=j+(1<<i);
			for(int k=j;k<l;k++){
				ll u=f[k],v=f[k|(1<<i)]*cur%MOD;
				f[k]=(u+v>MOD?u+v-MOD:u+v);
				f[k|(1<<i)]=(u>=v?u-v:u-v+MOD);
				cur=cur*step%MOD;
			}
		}
	}
	if(typ==-1){
		ll val=ksm(1<<l,MOD-2);
		for(int i=0;i<(1<<l);i++)
			f[i]=val*f[i]%MOD;
	}
	return;
}
poly operator *(poly f,poly g){
	int n=f.size()+g.size();
	poly h(n-1,0);
	int l=0;
	while(!g[l]) l++;
	for(int i=0;i<f.size();i++) if(f[i]) for(int j=l;j<g.size();j++)
		h[i+j]=(h[i+j]+f[i]*g[j])%MOD;
	return h;
}
}using namespace polynomial;
vector<int> son[MAXN];int dep[MAXN],siz[MAXN],n;
ll fac[MAXN],inf[MAXN],inv[MAXN],p[MAXN],res,sum;
ll C(ll x,ll y){return (x<0||x>y)?0:fac[y]*inf[x]%MOD*inf[y-x]%MOD;} 
poly dfs(int u){
	siz[u]=1;
	poly ans(dep[u]+2);
	for(int i=0;i<=dep[u]+1;i++) ans[i]=inf[i];
	for(int v:son[u]){
		dep[v]=dep[u]+1;
		ans=ans*dfs(v);
		siz[u]+=siz[v];
	}
	ans.resize(siz[u]+dep[u]+1);
	for(int i=0;i<siz[u];i++) ans[i]=0;
	return ans;
}
int main(){
	ios::sync_with_stdio(false);
	// freopen("Otomachi_Una.in","r",stdin);
	// freopen("Otomachi_Una.out","w",stdout);
	fac[0]=inf[0]=fac[1]=inf[1]=inv[1]=1;
	for(int i=2;i<MAXN;i++){
		fac[i]=fac[i-1]*i%MOD;inv[i]=(MOD-(MOD/i)*inv[MOD%i]%MOD);
		inf[i]=inf[i-1]*inv[i]%MOD; 
	}
	cin>>n;
	for(int i=2,f;i<=n;i++){
		cin>>f;son[f].push_back(i);
	}
	cout<<dfs(1)[n]*fac[n]%MOD<<'\n';
	return 0;
}