一、题目
洛谷的题目据说是转化过的，但是原来的题面太长我真不想看了。
二、解法
显然是两类元素问题，那么我们以白边为主，考虑原图会形成若干个置换环。
那么环内部是不能有任何白边的，然后我们把环当成点，不难发现问最后能形成多少个 ( t DAG)
补充：( t DAG) 计数是一类可以用 (dp) 解决的经典问题，基本问题是 (n) 个点有标号 ( t DAG) 计数。
设 (dp[i]) 表示 (i) 个点能构成的合法 ( t DAG)，每次我们枚举所有出度为 (0) 的点，就能找到子问题：
[dp[i]=sum_{j=1}^i dp[i-j]cdot {ichoose j}cdot 2^{j(i-j)} ]但是这样显然会算重，因为我们没有保证 (i-j) 中的点有可能度数为 (0)，考虑容斥去重，设钦定 (i) 个点出度为 (0) 的容斥系数为 (f_i)，那么考虑一个真实具有 (x) 个点的状态会被统计的次数：
[sum_{i=1}^x {xchoose i}cdot f_i=1 ]那么直接让 (f_i=(-1)^{i+1})，那么正确的转移应该是：
[dp[i]=sum_{j=1}^i(-1)^{j+1}cdot dp[i-j]cdot {ichoose j}cdot 2^{j(i-j)} ]
可以延用上面容斥的思路，然后设 (F[i][j]) 表示 (i) 个点想 (j) 个点连边的方案数（实际的点数）：
[dp[s]=sum_t dp[s-t]cdot (-1)^{|t|+1}cdot F[sz(s-t)][sz(t)]
]
其中 (sz(s)) 就表示集合 (s) 的置换环中包含的点数，现在的问题在于求解 (F[i][j])，考虑钦定 (j) 个点连向某个点的方案数是 (j!cdot (n-j-1)!)，也就是考虑这个点的排列中前面和后面的放置方法，根据乘法原理：
[F[i][1]=sum_{j=0}^i{ichoose j}cdot j!cdot (n-j-1)!
]
[F[i][j]=F[i][j-1]cdot F[i][1]
]
其实 (dp) 的那部分还可以优化，考虑我们只需要记录每个大小置换环的个数，这样状态数就小多了，转移的时候从每类置换环里钦定入度为 (0) 的即可。
三、总结
图计数问题和容斥很有关联，有正难则反法、待定容斥系数法。
#include <cstdio>
const int M = 45;
const int N = 1600;
const int MOD = 1e9+7;
#define int long long
int read()
{
int x=0,f=1;char c;
while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
return x*f;
}
int n,a[M],vis[M],s[M],b[M],fac[M];
int dp[N],sz[N],p[N],w[N],f[M][M],C[M][M];
void add(int &x,int y) {x=(x+y)%MOD;}
void init()
{
C[0][0]=fac[0]=1;
for(int i=1;i<=n;i++) fac[i]=fac[i-1]*i%MOD;
for(int i=1;i<=n;i++)
{
C[i][0]=1;
for(int j=1;j<=i;j++)
C[i][j]=(C[i-1][j-1]+C[i-1][j])%MOD;
}
for(int i=0;i<n;i++)
{
f[i][0]=1;
for(int j=0;j<=i;j++)
add(f[i][1],fac[j]*fac[n-j-1]%MOD*C[i][j]);
for(int j=2;j<=n;j++)
f[i][j]=f[i][j-1]*f[i][1]%MOD;
}
}
signed main()
{
n=read();init();
for(int i=1;i<=n;i++)
a[i]=read();
for(int i=1;i<=n;i++)
{
int x=i,sz=0;
while(!vis[x]) vis[x]=1,x=a[x],sz++;
b[sz]++;
}
s[0]=dp[0]=1;
for(int i=1;i<=n;i++) s[i]=s[i-1]*(b[i]+1);
for(int x=1;x<s[n];x++)
{
int t=1;p[1]=0;w[1]=MOD-1;
for(int i=1;i<=n;i++)
{
int tmp=t,v=(x/s[i-1])%(b[i]+1);
sz[x]+=i*v;
for(int j=1;j<=v;j++)
for(int h=1;h<=tmp;h++)
{
p[++t]=p[h]+j*s[i-1];
w[t]=w[h]*(j%2?MOD-1:1)%MOD*C[v][j]%MOD;
}
}
for(int i=2;i<=t;i++)
add(dp[x],dp[x-p[i]]*w[i]%MOD*
f[sz[x]-sz[p[i]]][sz[p[i]]]);
}
printf("%lld
",dp[s[n]-1]);
}