又是之前做过的题现在忘掉了。。。感谢whh巨佬提醒
首先要会读懂长长的题面。。初始时第(i)个人有物品(i)，之后你钦定了一个分配方案，如果一个集合的人在初始时能交换物品达到比你的方案更优的效果就不合法
或者直接看形式化题面更香
设(p_i)表示(i)的喜好排列，(a_i)表示你要分配给(i)的物品。
若(x)在(p_y)中位于(a_y)之前则(x)向(y)连一条红色有向边，同时(a_x)向(x)连蓝边
发现这个图合法的条件是不存在有红边的环。
原因我感性理解就是对于(x ightarrow y)蓝边表示把(x)给(y)不会比分配方案劣，而(x ightarrow y)红边表示把(x)给(y)更优了。
这样一个环上的人按边交换后得到了更优的方案，不合法。
这样我们要做的是把蓝边连好会得到几个环，再把这些环连成DAG，求所有图的贡献和。
设(f_i=i!(n-i+1)!)则一个图贡献为(prod f_{deg_i})其中(deg_i)表示(i)连入红边个数，也就是这个图对应排列方案数。
对于DAG计数类问题通常用状压DP，每次钦定一些入度(这里好像是出度。。)为(0)的点，容斥系数为((-1)^{|S|+1})
但是这题数据范围没法状压DP，考虑大小相同的环本质相同，这样极大地减小了状态数。
#include<bits/stdc++.h>
using namespace std;
#define fp(i,l,r) for(register int (i)=(l);(i)<=(r);++(i))
#define fd(i,l,r) for(register int (i)=(l);(i)>=(r);--(i))
#define fe(i,u) for(register int (i)=front[(u)];(i);(i)=e[(i)].next)
#define mem(a) memset((a),0,sizeof (a))
#define O(x) cerr<<#x<<':'<<x<<endl
#define int long long
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return x*f;
}
void wr(int x){
if(x<0)putchar('-'),x=-x;
if(x>=10)wr(x/10);
putchar('0'+x%10);
}
const int mod=1e9+7;
inline void tmod(int &x){x%=mod;}
inline int norm(int x){return (x%mod+mod)%mod;}
int n,sz[52],A[52],fac[52],C[52][52],tr[52][52],num[100010],dp[100010],prod[52],q[100010],cnt,g[100010];
bool vis[52];
main(){
n=read();fac[0]=1;fp(i,1,n)tmod(fac[i]=fac[i-1]*i);
fp(i,0,n){
C[i][0]=1;
fp(j,1,i)tmod(C[i][j]=C[i-1][j-1]+C[i-1][j]);
}
fp(i,0,n-1){
tr[0][i]=1;
fp(j,0,i)tmod(tr[1][i]+=C[i][j]*fac[j]%mod*fac[n-j-1]);
fp(j,2,n-i)tmod(tr[j][i]=tr[j-1][i]*tr[1][i]);
}
fp(i,1,n)A[i]=read();
fp(i,1,n)if(!vis[i]){
int x=i,d=0;
while(!vis[x])vis[x]=1,++d,x=A[x];
++sz[d];
}
prod[0]=dp[0]=1;g[0]=-1;
fp(i,1,n)prod[i]=prod[i-1]*(sz[i]+1);
fp(i,0,prod[n]-1){
q[cnt=1]=0;g[1]=-1;
fp(j,1,n){
int pre=cnt,t=(i/prod[j-1])%(sz[j]+1);
num[i]+=t*j;
fp(k,1,pre)fp(l,1,t){
q[++cnt]=q[k]+l*prod[j-1];
tmod(g[cnt]=g[k]*(l&1?-1:1)*C[t][l]);
}
}
fp(j,2,cnt)tmod(dp[i]+=dp[i-q[j]]*g[j]%mod*tr[num[q[j]]][num[i-q[j]]]);
}
printf("%lld
",norm(dp[prod[n]-1]));
return 0;
}