CF1466H - Finding satisfactory solutions
给出一个排列\(A\)，询问合法的数组\(p_i\)的个数，数组中每个元素\(p_i\)为一个排列。如果合法，当且仅当：不能找到一个排列\(B_i\)，使得存在一个集合\(S\)，满足：
- \(\forall i\in S,B_i\in S\)。
- \(\forall i\in S\)，在排列\(p_i\)中，不存在\(A_i\)在\(B_i\)左边出现。
- \(\exist i\in S\)，在排列\(p_i\)中，\(A_i\)在\(B_i\)右边出现。
\(n\le 40\)。
由于题面太长直接去看了洛谷的翻译，后来感觉不对劲才知道洛谷的翻译已经把算法的第一步说出来了……
可以把限制建成一个图：对于某个\(i\)，如果\(j\)在\(A_i\)前出现，则连黑边\((i,j)\)；连白边\((i,A_i)\)。
发现如果满足题目的限制，当且仅当：这个图中的所有环中不存在黑边。
理解：
由于\(\forall i\in S,B_i\in S\)，可以视作\(S\)为若干个环的集合。可以只考虑一个环。此时\((i,B_i)\)视作\(i\)的出边。
此时如果存在\(A_i\neq 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\)。
\(\binom{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\n",f[i]);
ll ans=f[p[n]-1];
ans=(ans+mo)%mo;
printf("%lld\n",ans);
return 0;
}
