CF1466H

给出一个排列(A)，询问合法的数组(p_i)的个数，数组中每个元素(p_i)为一个排列。如果合法，当且仅当：不能找到一个排列(B_i)，使得存在一个集合(S)，满足：

1. (forall iin S,B_iin S)。
2. (forall iin S)，在排列(p_i)中，不存在(A_i)在(B_i)左边出现。
3. (exist iin S)，在排列(p_i)中，(A_i)在(B_i)右边出现。

(nle 40)。

---

由于题面太长直接去看了洛谷的翻译，后来感觉不对劲才知道洛谷的翻译已经把算法的第一步说出来了……

可以把限制建成一个图：对于某个(i)，如果(j)在(A_i)前出现，则连黑边((i,j))；连白边((i,A_i))。

发现如果满足题目的限制，当且仅当：这个图中的所有环中不存在黑边。

理解：

由于(forall iin S,B_iin S)，可以视作(S)为若干个环的集合。可以只考虑一个环。此时((i,B_i))视作(i)的出边。

此时如果存在(A_i eq B_i)，即触犯了限制3，也就是环中存在黑边。

先连边((i,A_i))，由于(A_i)是排列所以形成的是若干个环。现在要往图中加入若干条边，使其不出现新的环。此时贡献为(prod out_i!(n-out_i-1)!)。

把环缩点，问题类似于DAG计数。然后状压DP。大小一样的环等价可以压在一起。设(f_S)表示(S)中的点组成的DAG的贡献。这里(S)用高维向量表示。状态不会很多，设(p(n))表示(n)划分成的可重集合中，不同的可重集合的个数最多是多少。(p(40)=1440)。

[f_S=sum_{Tsubseteq S} f(Ssetminus T)inom{S}{T}(-1)^{|T|+1}calc(sum(T),sum(Ssetminus T))\ calc(p,q)=(sum_{i} inom{q}{i} i!(n-1-i)!)^p ]

(inom{S}{T})就是各个维度的组合数乘积。

直接这样算大概是(O(p^2(n)n))的。

这个东西可以辅助转移。考虑到(T)转移到(T+Delta)，可以给(Delta)一位位分别计算。另外设个(DP)存下(T+Delta)，当前处理到第几位，以及(sum_T)是多少（用于计算(calc)产生的贡献）。设(g_{S,i,j})，具体见程序。时间(O(p(n)n^3))。

---

using namespace std;
#include <bits/stdc++.h>
#define N 45
#define ll long long
#define mo 1000000007
ll qpow(ll x,ll y=mo-2){
	ll r=1;
	for (;y;y>>=1,x=x*x%mo)
		if (y&1)
			r=r*x%mo;
	return r;
}
int n;
int a[N];
int c[N],p[N];
ll fac[N*2],ifac[N*2];
void initC(int n){
	fac[0]=1;
	for (int i=1;i<=n;++i)
		fac[i]=fac[i-1]*i%mo;
	ifac[n]=qpow(fac[n]);
	for (int i=n-1;i>=0;--i)
		ifac[i]=ifac[i+1]*(i+1)%mo;
}
ll C(int m,int n){return fac[m]*ifac[n]%mo*ifac[m-n]%mo;}
void init(){
	scanf("%d",&n);
	initC(n*2);
	for (int i=1;i<=n;++i)
		scanf("%d",&a[i]);
	static int bz[N];
	for (int i=1;i<=n;++i)
		if (!bz[i]){
			int t=0;
			for (int x=i;!bz[x];x=a[x])
				bz[x]=1,t++;
			c[t]++;
		}
	p[0]=1;
	for (int i=1;i<=n;++i)
		p[i]=p[i-1]*(c[i]+1);
}
ll q[N];
ll f[1440];
ll g[1440][N][N];
void doit(){
	for (int i=0;i<=n;++i){
		ll s=0;
		for (int j=0;j<=i;++j)
			(s+=C(i,j)*fac[j]%mo*fac[n-1-j])%=mo;
		q[i]=s;
	}
	for (int i=0;i<p[n];++i){
		for (int t=0;t<=n;++t)
			for (int j=1;j<=n;++j){
				int w=i%p[j]/p[j-1];
				ll qj=qpow(q[t],j),tmp=1;
				for (int k=w;k<=c[j];++k){
					(g[i+(k-w)*p[j-1]][j][t]+=g[i][j-1][t]*ifac[k-w]%mo*tmp)%=mo;
					tmp=tmp*(-1)*qj%mo;
				}
			}
		ll S=1,s=0;
		if (i==0)
			f[0]=1;
		else{
			for (int j=1;j<=n;++j){
				int w=i%p[j]/p[j-1];
				S=S*fac[w]%mo;
				s+=w*j;
			}
			ll tmp=0;
			for (int j=0;j<=s;++j)
				(tmp+=g[i][n][j])%=mo;
			f[i]=tmp*S%mo;
		}
//		(g[i][0][s]+=f[i]*qpow(S)*(-1))%=mo;
		ll v=f[i]*qpow(S)*(-1)%mo;
		for (int j=1;j<=n;++j){
			int w=i%p[j]/p[j-1];
			ll qj=qpow(q[s],j),tmp=(-1)*qj;
			for (int k=w+1;k<=c[j];++k){
				(g[i+(k-w)*p[j-1]][j][s]+=v*ifac[k-w]%mo*tmp)%=mo;
				tmp=tmp*(-1)*qj%mo;
			}
		}
	}
}
int main(){
	freopen("in.txt","r",stdin);
	init();
	doit();
//	for (int i=0;i<p[n];++i)
//		printf("%lld
",f[i]);
	ll ans=f[p[n]-1];
	ans=(ans+mo)%mo;
	printf("%lld
",ans);
	return 0;
}