主要思路
这是道好题，运用了图论建模的方法，可谓是十分精妙了。
对于每一个\(i\)与\(arr[i]\)都连一条有向边，比如说数列\(2\ 1\ 4\ 3\)：
可以通过 DFS 染色的方式找出\(k\)个环，并且将每个环的长度记为\(l_i\)。我们可以得出一个引理：
在图中，对于长度为\(n\)的环，将其变为自环的最小步骤为\(n-1\)。
就是这样的性质，我们可以记\(F[i]\)为当环长为\(i\)时变为自环的方案数，那么我们可以基于分裂的思想，设\(T(x,y)\)为当环长为\(x+y\)时分裂为环长分别为\(x,y\)时的方案数，显然：
\[T(x,y) = \begin{cases} \frac{n}{2}, x = y \\ n, x \neq y \end{cases}\]
所以，我们可以照例推出：
\[ F[i] = \sum_{x+y=i} T(x,y)*F[x]*F[y]*\frac{(n-2)!}{(x-1)!(y-1)!} \]
我来解释一下这个式子的意义：首先，枚举\(x\)和\(y\)的情况，然后乘上这两个本身变成自环的方案数，也就是\(F[x]*F[y]\)，之后根据多重集排列公式和上面那个引理，两者步数分别为\(x-1,y-1\)，所以也不难得出上面那一坨了。
答案为\(\prod_{i = 1}^k {F_{l_i}} * \frac{(n-k)!}{\prod_{i = 1}^{k}(l_i-1)}\)，时间复杂度为\(O(n^2)\)，然而通过人类智慧等方法，发现\(F[n]=n^{n-2}\)，所以通过快速幂降为\(O(n \log n)\)。
代码
// CH3602.cpp
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int MAX_N = 1e5 + 200, mod = 1e9 + 9;
int T, n, arr[MAX_N], head[MAX_N], current, k, li[MAX_N], vis[MAX_N];
ll level[MAX_N], inv[MAX_N];
struct edge
{
int to, nxt;
} edges[MAX_N << 1];
void addpath(int src, int dst)
{
edges[current].to = dst, edges[current].nxt = head[src];
head[src] = current++;
}
int dfs(int u, int fa)
{
vis[u] = true;
for (int i = head[u]; i != -1; i = edges[i].nxt)
if (edges[i].to != fa)
{
if (vis[edges[i].to])
return 1;
return dfs(edges[i].to, u) + 1;
}
return 1;
}
ll quick_power(ll bas, ll tim)
{
ll ans = 1;
while (tim)
{
if (tim & 1)
ans = ans * bas % mod;
bas = bas * bas % mod;
tim >>= 1;
}
return ans;
}
ll Fn(int n) { return (n == 1) ? 1 : quick_power(n, n - 2); }
int main()
{
scanf("%d", &T);
level[0] = inv[0] = 1;
for (int i = 1; i < MAX_N; i++)
level[i] = level[i - 1] * i % mod;
inv[MAX_N - 1] = quick_power(level[MAX_N - 1], mod - 2);
for (int i = MAX_N - 1; i >= 2; i--)
inv[i - 1] = inv[i] * i % mod;
while (T--)
{
memset(head, -1, sizeof(head)), current = 0;
memset(li, 0, sizeof(li)), k = 0;
memset(vis, 0, sizeof(vis));
scanf("%d", &n);
for (int i = 1; i <= n; i++)
scanf("%d", &arr[i]), addpath(arr[i], i), addpath(i, arr[i]);
for (int i = 1; i <= n; i++)
if (!vis[i])
li[++k] = dfs(i, 0);
ll ans = 1;
for (int i = 1; i <= k; i++)
ans = ans * Fn(li[i]) % mod;
ans = (ans * level[n - k] % mod);
for (int i = 1; i <= k; i++)
ans = ans * inv[li[i] - 1] % mod;
printf("%lld\n", ans);
continue;
}
return 0;
}