提交记录

大佬们帮忙看看本蒟蒻的代码哪里错了？（一脸茫然）

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 985661441;
const int maxn1 = 3000005;
const int maxn2 = 5005;
ll qpzc(ll a, ll b) {
	if (b < 0) return qpzc(qpzc(a, mod - 2), -b);
	ll ans = 1;
	while (b) {
		if (b & 1) (ans *= a) %= mod;
		(a *= a) %= mod;
		b >>= 1;
	}
	return ans;
}
ll fac[maxn1], inv[maxn1];
void orz() {
	fac[0] = 1;
	for (ll i = 1; i <= maxn1 - 5; i++) fac[i] = fac[i - 1] * i % mod;
	inv[maxn1 - 5] = qpzc(fac[maxn1 - 5], mod - 2);
	for (ll i = maxn1 - 6; i >= 0; i--) inv[i] = inv[i + 1] * (i + 1) % mod;
}
ll C(ll i, ll j) {
	if (i < 0 || j < 0 || i < j) return 0;
	return fac[i] * inv[j] % mod * inv[i - j] % mod;
}
vector<ll> vc[maxn2];
ll dp[maxn2][maxn2][3];
ll siz[5005], tmp[5005][3], tot[5005];
void dfs(ll now) {
	siz[now] = 1;
	dp[now][0][0] = 1;
	for (auto v : vc[now]) {
		dfs(v);
		for (ll i = 0; i <= siz[now] + siz[v]; i++) tmp[i][0] = tmp[i][1] = tmp[i][2] = 0;
		for (ll i = 0; i < siz[now]; i++) {
			for (ll j = 0; j < siz[v]; j++) {
				ll g = (1ll * dp[v][j][0] + dp[v][j][1] + dp[v][j][2]) % mod;
				(tmp[i + j + 1][0] += 1ll * dp[now][i][0] * g) %= mod;
				(tmp[i + j + 1][1] += 1ll * dp[now][i][1] * g) %= mod;
				(tmp[i + j + 1][2] += 1ll * dp[now][i][2] * g) %= mod;
				(tmp[i + j][1] += 1ll * dp[now][i][0] * (1ll * dp[v][j][0] + dp[v][j][1])) %= mod;
				(tmp[i + j][2] += 1ll * dp[now][i][1] * (1ll * dp[v][j][0] + dp[v][j][1])) %= mod;
			}
		}
		siz[now] += siz[v];
		for (ll i = 0; i < siz[now] + siz[v]; i++) dp[now][i][0] = tmp[i][0], dp[now][i][1] = tmp[i][1], dp[now][i][2] = tmp[i][2];
	}
}
int main() {
	orz();
	ll n, ans = 0;
	cin >> n;
	for (ll i = 2; i <= n; i++) {
		int f;
		cin >> f;
		vc[f].push_back(i);
	}
	dfs(1);
	for (int i = 0; i < n; i++) {
		tot[i] = (1ll * dp[1][i][0] + dp[1][i][1] + dp[1][i][2]) % mod * inv[n - 1] % mod * fac[n - 1 - i] % mod * fac[i] % mod;
		for (ll j = 0; j < i; j++) (tot[i] += mod - tot[j]) %= mod;
		(ans += tot[i] * i % mod) %= mod;
	}
	cout << ans;
	return 0;
}