この問題は、「すべての要素が特定の条件を満たさない」ような場合の数え上げ問題であり、直接的に解くか、または「少なくとも1つの要素が条件を満たす」場合に包除原理を適用する方法があります。ここでは後者のアプローチを取ります。
(S(i)) を、「(i) 個の (k) について (a_k) が (k) に到達可能である」という条件を満たす順列の数と定義します。包除原理を使うと、以下のように表現できます。
[ans = \sum_{i=0}^n (-1)^i S(i) ]
(S(i)) を求めるために、動的計画法 (DP) を使用します。(f_{u, i}) を、「(u) の部分木内に (i) 個の (k) があり、それらの (a_k) が (k) に到達可能である」ような場合の数と定義します。ここで、(siz_u) は (u) の部分木（(u) 含む）のサイズです。
まず、以下の2つのケースを考えます：
- (u) が (u) の部分木（(u) 除外）内の任意の (k) の (a_k) である場合。
- (u) が (u) の部分木（(u) 除外）内の任意の (k) の (a_k) でない場合。
ケース1: まず、ケース2の部分を求めます。その後、(siz_u) から 1 まで降順にスキャンし、以下の更新を行います：
[\forall i \in [1, siz_u], f_{u, i} \gets f_{u, i} + f_{u, i-1} \times ((siz_u - 1) - (i - 1)) ]
ここで、(siz_u - 1) は (u = a_k) となるときの選択可能な (k) の数、(i - 1) は (u) の子供たちが既に選んだ (k) の数です。
ケース2: 各子供 (v) に対して、その貢献度を計算します。最初の考え方は複雑ですが、効率的な方法を見つけることができます。(u) の子供を (v_1, \ldots, v_s) として、(sf_{p-1, i}) を前 (p-1) 個の子供が (i) 個の (a_k) が (k) に到達可能である場合の数と定義します。すると、次の式で更新できます：
[ sf_{p, i} = \sum_{j=0}^i f_{v, j} \times sf_{p-1, i-j} ]
最終的には (f_{u, i} = sf_{s, i}) となります。具体的な実装は以下の通りです：
int f[N][N], cf[N], siz[N];
void fdp(int u) {
f[u][0] = siz[u] = 1;
for (int v : vc[u]) {
fdp(v);
memcpy(cf, f[u], (siz[u] + 1) << 2);
memset(f[u], 0, (siz[u] + 1) << 2);
for (int i = 0; i <= siz[u]; ++i) {
for (int j = 0; j <= siz[v]; ++j) {
(f[u][i + j] += 1LL * cf[i] * f[v][j] % mod) %= mod;
}
}
siz[u] += siz[v];
}
for (int i = siz[u]; i > 0; --i) {
(f[u][i] += 1LL * f[u][i - 1] * (siz[u] - i) % mod) %= mod;
}
}
この方法は、各要素が前方のすべての要素とのみ関連を持つようにすることで重複を防ぎ、全体の計算量を (\Theta(n^2)) に抑えることができます。
最後に、残りの (n - i) 個の位置は自由に並べられるため、(S(i) = f_{1, i} \times (n - i)!) となります。
以上のコードの計算量は、以下のようになります：
[\sum siz_u \sum_{v \in son(u)} siz_v = \sum (siz_u)^2 = \Theta(n^3)]
これを改善するために、(v_p) に対する (sf_{p, i-1}) で (i > \sum_{p'=1}^{p-1} siz(v_{p'})) の部分は無視できるという点を利用します。siz[u] を動的に更新し、(i) を 0 から (\sum_{p'=1}^{p-1} siz(v_{p'})) まで枚挙します：
int f[N][N], cf[N], siz[N];
void fdp(int u) {
f[u][0] = siz[u] = 1;
for (int v : vc[u]) {
fdp(v);
memcpy(cf, f[u], (siz[u] + 1) << 2);
memset(f[u], 0, (siz[u] + 1) << 2);
for (int i = 0; i <= siz[u]; ++i) {
for (int j = 0; j <= siz[v]; ++j) {
(f[u][i + j] += 1LL * cf[i] * f[v][j] % mod) %= mod;
}
}
siz[u] += siz[v];
}
for (int i = siz[u]; i > 0; --i) {
(f[u][i] += 1LL * f[u][i - 1] * (siz[u] - i) % mod) %= mod;
}
}
この改善により、計算量は (\Theta(n^2)) になります。
包除原理を使用する理由は、問題の制約（(a) が ({1, \ldots, n}) の順列であること）を満たすためです。DP の状態設計では、(a_k) が部分木内に制限されますが、部分木は互いに交わらないため、重複が発生しません。逆に、条件を満たさない (a_k) は根へのパス上にあるため、部分木での DP が容易になります。