题意：给定一个 1∼n1\sim n 的排列 p1,p2,…,pnp_1,p_2,…,p_n，可进行若干次操作，每次选择两个整数 x,yx,y，交换 px,pyp_x,p_y。

设把 p1,p2,…,pnp_1,p_2,…,p_n 变成单调递增的排列 1,2,…,n1,2,…,n 至少需要 mm 次交换。

求有多少种操作方法可以只用 mm 次交换达到上述目标。

因为结果可能很大，你只需要输出结果对 109+910^9+9 取模之后的值。

例如排列 2,3,12,3,1 至少需要 22 次交换才能变为 1,2,31,2,3。操作方法共有 33 种，分别是：

方法一：先交换数字 2,32,3，变成 3,2,13,2,1，再交换数字 3,13,1，变成 1,2,31,2,3。
方法二：先交换数字 2,12,1，变成 1,3,21,3,2，再交换数字 3,23,2，变成 1,2,31,2,3。
方法三：先交换数字 3,13,1，变成 2,1,32,1,3，再交换数字 2,12,1，变成 1,2,31,2,3。

对于每个排列 p1,p2,…,pnp_1,p_2,…,p_n，我们从每个 ii 向 pip_i 连一条边，可以得到 nn 个点 nn 条边的图，由若干个环构成。

最后目标排列为 1,2,⋅⋅⋅,n1,2,···,n ，由 nn 个自环构成。

引理：

把一个长度为 nn 的环变成 nn 个自环，最少需要有 n−1n-1 次操作。

设 FnF_n 表示用最少的步数把一个长度为 nn 的环变成 nn 个自环的操作方法数。

在操作中，我们可以把该环拆成长度为 xx 和 yy 的两个环，其中 x+y=nx+y=n。设 T(x,y)T(x,y) 表示有多少种交换方法可以把长度为 nn 的环变成长度位 x,yx,y 的两个环，容易发现：

T(x,y)={n2n 是偶数且 x=ynn 是奇数或 x≠y T(x,y)= \begin{cases} \frac{n}{2}& n \text{ 是偶数且 } x=y \\ n& n \text{ 是奇数或 } x \ne y \\ \end{cases}

另外，二者各自变为自环的方法数为 Fx,FyF_x,F_y，步数为 x−1x-1 和 y−1y-1。

可得：

Fn=∑x+y=nT(x,y)×Fx×Fy×(n−2)!(x−1)!(y−1)! F_n=\sum_{x+y=n}T(x,y)\times F_x\times F_y\times\frac{(n-2)!}{(x-1)!(y-1)!}

如果最初的排列构成的图，有长度为 l1,l2,⋅⋅⋅lkl_1,l_2,···l_k 的 kk 个环组成，所以最终答案为：

Fl1×Fl2×⋅⋅⋅×Flk×(n−k)!(l1−1)!(l2−1)!⋅⋅⋅(lk−1)! F_{l_1}\times F_{l_2}\times ···\times F_{l_k}\times \frac{(n-k)!}{(l_1-1)!(l_2-1)!···(l_k-1)!}

于是bfs找环解决，事实上，我们可以发现，Fn=nn−2F_n=n^{n-2}，可以进一步优化。

code:

#include<cstdio>
#include<cmath>
#include<cstring>
#include<cstdlib>
#include<algorithm>
using namespace std;
const int mod=1e9+9;
int t,n,tot;
int head[401010],cnt;
long long fac[101010],inv[101010],ans;
int vise[401010],visv[101010],len[101010];
struct node
{
	int next;
	int to;
}e[401010];
void add(int from,int to)
{
	e[++cnt].next=head[from];
	e[cnt].to=to;
	head[from]=cnt;
}
void dfs(int u,int p,int r)
{
	if(!visv[u])
	{
		len[p]++;
		visv[u]=1;
	}
	for(int i=head[u];i;i=e[i].next)
	{
		if(vise[i]==r)
			continue;
		vise[i]=vise[i^1]=r;
		dfs(e[i].to,p,r);
	}
}
long long power(long long a,long long b,long long p)
{
	if(b<=0)
		return 1;
	long long ans=1%p;
	for(;b;b>>=1)
	{
		if(b&1)
			ans=ans*a%p;
		a=a*a%p;
	}
	return ans%p;
}
int main()
{
	fac[0]=inv[0]=1;
	for(int i=1;i<=100000;i++)
	{
		fac[i]=(fac[i-1]%mod*i%mod)%mod;
		inv[i]=power(fac[i],mod-2,mod)%mod;
	}
	scanf("%d",&t);
	for(int _=1;_<=t;_++)
	{
		memset(e,0,sizeof(e));
		memset(vise,0,sizeof(vise));
		memset(visv,0,sizeof(visv));
		memset(len,0,sizeof(len));
		memset(head,0,sizeof(head));
		ans=0;
		cnt=1;
		tot=0;
		scanf("%d",&n);
		for(int i=1;i<=n;i++)
		{
			int x;
			scanf("%d",&x);
			add(i,x);
			add(x,i);
		}
		for(int i=1;i<=n;i++)
			if(!visv[i])
				dfs(i,++tot,_);
		ans=fac[n-tot];
		for(int i=1;i<=tot;i++)
			ans=(ans%mod*power(len[i],len[i]-2,mod)%mod*inv[len[i]-1]%mod)%mod;
		printf("%lld\n",ans);
	}
    return 0;
}