给出一个排列(A)，询问合法的数组(p_i)的个数，数组中每个元素(p_i)为一个排列。如果合法，当且仅当：不能找到一个排列(B_i)，使得存在一个集合(S)，满足：
- (forall iin S,B_iin S)。
- (forall iin S)，在排列(p_i)中，不存在(A_i)在(B_i)左边出现。
- (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;
}