考虑如果确定了 \(a\)，怎么去算 \(b\) 的方案数。
对于一次赋值 \(b_{p_i}\leftarrow p_i\)，考虑钦定其为最后一次对 \(b_{a_{p_i}}\) 的修改，那么相当于把 \(a_{p_i}\) 锁定。
于是，如果当前 \(p_i=x\)，那么我们可以决定是否把 \(x\) 锁定，以及是否把 \(a_x\) 锁定。
考虑如果在 \(x\) 把 \(a_x\) 锁定，那么将 \(x\) 向 \(a_x\) 连一条边。显然新图由若干条链和单点组成，而且我们可以通过新图还原出唯一一个 \(b\)。
但是这样的链划分不一定全部合法。对于一个单点 \(x\)，需要 \(a_x\) 不为单点，且 \(a_x\) 不为一条链的起点。这个限制也是充分的。
这个做法可以对每个 \(a\) 求答案，但是并不好做对所有 \(a\) 求和。考虑反过来，对于一个已知的新图，算合法的 \(a\) 个数。
设新图中有 \(i\) 个单点，\(j\) 条链，那么首先组合出这个图的方案数是 \(\dbinom{n}{i}\dfrac{(n-i)!}{j!}\dbinom{n-i-j-1}{j-1}\)。然后对于单点，其出边必须是链上点，且不为起点，方案数 \((n-i-j)^i\)。对于链，可以任意连，方案数 \((n-1)^j\)。乘起来即可。
#include <bits/stdc++.h>
using namespace std;
using ll = long long;int mod;
void Add(int &x, ll y) {x = (x + y) % mod;
}
int Pow(int x, int y) {int b = x, r = 1;for(; y; b = (ll)b * b % mod, y /= 2) {if(y & 1) r = (ll)r * b % mod;}return r;
}const int kN = 5005;
int n;
int mul[kN], imul[kN], pw[kN][kN];void Init(int N = kN - 2) {mul[0] = 1;for(int i = 1; i <= N; i++) {mul[i] = (ll)mul[i - 1] * i % mod;}imul[N] = Pow(mul[N], mod - 2);for(int i = N - 1; ~i; i--) {imul[i] = (ll)imul[i + 1] * (i + 1) % mod;}for(int i = 0; i <= n; i++) {pw[i][0] = 1;for(int j = 1; j <= n; j++) {pw[i][j] = (ll)pw[i][j - 1] * i % mod;}}
}
int C(int n, int m) {if((n < m) || (m < 0)) {return 0;}return (ll)mul[n] * imul[m] % mod * imul[n - m] % mod;
}int main() {// freopen("1.in", "r", stdin);// freopen("1.out", "w", stdout);ios::sync_with_stdio(0), cin.tie(0);cin >> n >> mod;Init();int res = 0;for(int i = 0; i <= n; i++) {for(int j = 1; j * 2 + i <= n; j++) {int val = C(n, i);val = (ll)val * mul[n - i] % mod * imul[j] % mod * C(n - i - j - 1, j - 1) % mod;val = (ll)val * pw[n - i - j][i] % mod;val = (ll)val * pw[n - 1][j] % mod;Add(res, val);}}cout << res << "\n";return 0;
}