题意简述
有 (n) 个 (1) 到 (n) 的排列，第 (i) 个排列为 (P_i)。
定义一个排列 (A) 是最优的，当且仅当不存在任意一个排列 (B) 满足存在一个下标集合 (S) 符合以下所有条件：
（1）(forall iin S)，(B_iin S)。
（2）(forall iin S)，在 (P_i) 中 (A_i) 不能出现在 (B_i) 的左边。
（3）(exists iin S)，在 (P_i) 中 (B_i) 出现在 (A_i) 的左边。
给定一个排列 (A)，求有多少种 (P_{1dots n}) 使得排列 (A) 是最优的，对 (10^9+7) 取模。
(1le nle 40)。
Solution
显然在 (P_i) 中只有出现在 (A_i) 及之前的数是有用的。
考虑建一张有向图：如果在 (P_i) 中 (j) 出现的位置严格在 (A_i) 出现的位置左边，则连白边 ((i,j))，然后对于所有的 (i) 连黑边 ((i,A_i))。
对应到排列 (A) 不是最优的条件，相当于可以选出若干个不相交的环（符合了条件 (1)），满足这些环的边不全为黑（符合了条件 (3)）。
到这里可以发现其实没有必要选若干个不相交的环，只要能选一个不全为黑边的环出来，(A) 就不是最优的了。
考虑一个强连通分量，不难发现只要分量内有白边，就一定能选出一个不全为黑边的环出来。
于是我们得到 (A) 是最优的条件：图的每个强连通分量内都只有黑边。
由于所有的黑边 ((i,A_i)) 已经给出，并且所有黑边构成了若干个不相交的环，故把这些环当成一个大点，要求的就是把这些大点连成一个 DAG 的方案数。
同时由于 (P_i) 是排列，故如果 (i) 点出发的白边有 (d) 条，则会带上 (d!(n-d-1)!) 的乘积贡献（把 (P_i) 中 (A_i) 之前和之后的数安排顺序）。
考虑暴力状压 DP：(f_S) 表示把集合 (S) 内的大点连成 DAG 的方案数。
转移考虑枚举入度为 (0) 的点集 (T)，加上点集 (T) 连向点集 (S-T) 的方案数乘上 (f_{S-T})。
这样不能保证 (S-T) 内没有入度为 (0) 的点，故还要套一个容斥：设上面算出的结果为 (g_T)，则真正入度为 (0) 的点集恰好为 (T) 的方案数应为 (sum_{Tsubseteq Usubseteq S}(-1)^{|U|-|T|}g_U)。
于是：
这样就可以枚举子集进行暴力 DP 了，复杂度 (O(2^n))。
而本题中我们可以把 (S) 记录成每种大小的大点分别出现的次数，枚举子集可以改成枚举 (U) 中每种大小的大点出现了多少次后乘一个组合数，这样处理之后复杂度可以通过。
Code
#include <bits/stdc++.h>
#define mp std::make_pair
typedef std::pair<int, int> pii;
typedef std::vector<pii> vpii;
template <class T>
inline void read(T &res)
{
res = 0; bool bo = 0; char c;
while (((c = getchar()) < '0' || c > '9') && c != '-');
if (c == '-') bo = 1; else res = c - 48;
while ((c = getchar()) >= '0' && c <= '9')
res = (res << 3) + (res << 1) + (c - 48);
if (bo) res = ~res + 1;
}
const int N = 45, M = 3e5 + 5, EI = 1e9 + 7;
int n, p[N], ToT, son[M][N], Root, f[M], C[N][N], g[N], h[N][N], fac[N];
vpii nums[M];
bool vis[N];
int dfs(int cur, vpii now)
{
int u; nums[u = ++ToT] = now; vpii nxt = now;
int lst = now.empty() ? 0 : now[now.size() - 1].first;
for (int i = lst; i <= n - cur; i++) if (i)
{
if (i == lst) nxt[nxt.size() - 1].second++;
else nxt.push_back(mp(i, 1));
son[u][i] = dfs(cur + i, nxt);
nxt = now;
}
return u;
}
void dfs_trans(int u, int dep, int v, int in, int out,
int now, vpii &trans)
{
if (dep == nums[u].size())
{
if (u == v) return;
trans.push_back(mp(1ll * now * h[out][in] % EI, v));
return;
}
int x = nums[u][dep].second, y = nums[u][dep].first;
for (int i = 0; i <= x; i++)
dfs_trans(u, dep + 1, v, in, out, 1ll * now * C[x][i] % EI, trans),
v = son[v][y], in += y, out -= y, now = EI - now;
}
int dp(int u)
{
if (f[u] != -1) return f[u];
if (u == Root) return f[u] = 1;
vpii trans; int sze = 0, cnt = 0; f[u] = 0;
for (int i = 0; i < nums[u].size(); i++)
sze += nums[u][i].first * nums[u][i].second,
cnt += nums[u][i].second;
dfs_trans(u, 0, Root, 0, sze, cnt & 1 ? 1 : EI - 1, trans);
for (int i = 0; i < trans.size(); i++)
f[u] = (1ll * trans[i].first * dp(trans[i].second) + f[u]) % EI;
return f[u];
}
int main()
{
read(n); Root = dfs(0, vpii(0));
for (int i = 1; i <= n; i++) read(p[i]);
std::vector<int> cir;
for (int i = 1; i <= n; i++) if (!vis[i])
{
int c = 0;
for (int j = i; !vis[j]; j = p[j]) c++, vis[j] = 1;
cir.push_back(c);
}
int u = 1;
std::sort(cir.begin(), cir.end());
for (int i = 0; i < cir.size(); i++) u = son[u][cir[i]];
for (int i = 0; i <= n; i++) C[i][0] = 1;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= i; j++)
C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % EI;
memset(f, -1, sizeof(f)); fac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = 1ll * fac[i - 1] * i % EI;
for (int i = 0; i < n; i++)
{
for (int j = 0; j <= i; j++)
g[i] = (1ll * fac[j] * fac[n - j - 1] % EI * C[i][j] + g[i]) % EI;
h[0][i] = 1;
for (int j = 1; j <= n; j++)
h[j][i] = 1ll * h[j - 1][i] * g[i] % EI;
}
return std::cout << dp(u) << std::endl, 0;
}