C.拓扑

考虑dp，如果以子树内的状态来更新的话，实际上就是考虑子树内有多少个点的拓扑序在当前树根之前，但是会发现这样实际上会影响根在拓扑序中的位置，这是不好维护的，所以考虑定义dp状态为:$f_{i,j}$表示$i$以及$i$子树外的节点的相对顺序都处理好了，并且$p_{i}=i$的方案数，这样定义无论$i$子树内部是怎么放，$p_{i}=i$都是恒成立的。

考虑转移，也就是$f_{p,j}\Rightarrow f_{u,i}$，其中$u$是$p$的儿子，此时我们需要枚举有多少个$u$的兄弟子树的节点排在$i$（也就是$u$在拓扑序中的位置）之前，这个数量记为$k$，同时需要保证$p$的兄弟子树的节点中有$x$个节点也在$i$之前，所以有$x+k=i-j-1$，注意状态定义里强调的是相对位置而不是绝对位置，所以说这一部分就是把$k$个节点插入到这$x+1$个缝隙当中，用挡板法很容易得到方案数为$C_{x+k}^{k}=C_{i-j-1}^{k}$，右边同理，但是要注意，我们没有处理$u$的子树中（不包含$u$）的点，所以要预留$sz_{u}-1$个位置，所以这一部分对应的方案数为$C_{n-i-sz_{u}+1}^{sz_{p}-1-sz_{u}-k}$，至此，我们可以得到转移式：

$f_{u,i}=\sum_{j=1}^{i-1}f_{p,j}\times \sum_{k=0}^{i-j-1} C_{i-j-1}^{k}\times C_{n-i-sz_{u}+1}^{sz_{p}-1-sz_{u}-k} \times tmp$ 其中$tmp$是指$u$的兄弟子树中的点的拓扑序的方案数，至于怎么算放到后面介绍

这个转移式的复杂度过高，我们发现后面的组合数相乘，可以用范德蒙德卷积优化，化简后得到：

$f_{u,i}=\sum_{j=1}^{i-1}f_{p,j}\times C_{n-sz_{u}-j}^{sz_{p}-1-sz_{u}} \times tmp$，这样我们就可以用前缀和优化这个dp了，时间复杂度$O(N^2)$

最后讲一下tmp怎么得到，我们定义$dp_i$表示只考虑$i$ 子树内部的点所构成的拓扑序的方案数，那么转移即为$dp_u=(\prod_{v \in son(u)} dp_v) \times \frac{(sz_u-1)!}{\prod_{v \in son(u)} sz_v!}$

那么$tmp$就是消除$u$的子树对$dp_p$的贡献，很容易得到$tmp=\frac{dp_p}{dp_u\times C_{sz_p-1}^{sz_u}}$

最后在输出答案的时候，别忘了把$i$子树内部的点插入到最后，具体见代码。

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define ll long long
#define PII pair<int,int>
#define lson 2*p
#define rson 2*p+1
#define x first
#define y second
#define LL long long
#define lowbit(x) x&-x
#define endl "\n"
const int N=5e3+10;
const int mod=998244353;
vector<int> g[N];
int n,f[N][N],p[N];
int fact[N],infact[N],sz[N];
int dp[N];
int qpow(int a,int b) {
    int ans=1;
    while(b) {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
void init() {
    fact[0]=1;
    for(int i=1;i<N;i++) fact[i]=fact[i-1]*i%mod;
    infact[N-1]=qpow(fact[N-1],mod-2);
    for(int i=N-2;i>=0;i--) infact[i]=infact[i+1]*(i+1)%mod;
}
int c(int n,int m) {
    if(n<m) return 0;
    return fact[n]*infact[m]%mod*infact[n-m]%mod;
}
void dfs(int x,int fa) {
    sz[x]=dp[x]=1;
    int res=1;
    for(int y:g[x]) {
        if(y==fa) continue;
        dfs(y,x);
        sz[x]+=sz[y];
        dp[x]=dp[x]*dp[y]%mod;
        res=res*infact[sz[y]]%mod;
    }
    dp[x]=dp[x]*fact[sz[x]-1]%mod*res%mod;
}
void dfs1(int p,int fa) {
    for(int u:g[p]) {
        if(u==fa) continue;
        int res=dp[p]*qpow(dp[u]*c(sz[p]-1,sz[u])%mod,mod-2)%mod;
        int sum=0;
        for(int i=1;i<=n;i++) {
            sum+=f[p][i-1]*c(n-i+1-sz[u],sz[p]-sz[u]-1)%mod*res%mod;
            sum%=mod;
            f[u][i]+=sum%mod;
            f[u][i]%=mod;
        }
        dfs1(u,p);
    }
}
void solve() {
    cin>>n;
    for(int i=2;i<=n;i++) {
        cin>>p[i];
        g[p[i]].push_back(i);
        g[i].push_back(p[i]);
    }
    dfs(1,0);
    f[1][1]=1;
    dfs1(1,0);
    for(int i=1;i<=n;i++) {
        int ans=f[i][i]*dp[i]%mod*c(n-i,sz[i]-1)%mod;
        cout<<ans<<" ";
    }
    cout<<endl;
}  
signed main() {
    int T=1;
    ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
    init();
    // cin>>T;
    while(T--) {    
        solve();
    }
}