非常 adhoc 的题，但是值得一练的好题！
一眼下去，我们会发现这个条件真的太过于抽象，根本无法想象。
注意到题目给了我们一个关键信息：一个排列组 \(\{b_1,b_2,...,b_n\}\) 对应唯一的好的分配方案。
考虑建立图论模型：每个人的编号向最喜欢的物品编号连边，形成一棵内向基环树森林。
对于基环树森林中的一个环，我们发现一定对应的 \(a\) 中的这个环。如果这个环没有在 \(a\) 中出现，那么环上每个点调整成这个环后会更优。
所以我们可以从森林中删掉这个环，然后继续一直这样做下去，最终会形成若干个环，就是 \(a\)。
考虑删去一个环后，原来连边指向这个环上的点的所有点都会重新调整连边，即连向 \(b_i\) 中下一个没有被删掉的点。
因此对于 \(i\)，设 \(a_i\) 在 \(b_i\) 中的出现位置为 \(p\)，那么 \(b_{i,1...p-1}\) 都是先前被删去的点。
这样对 \(a\) 中的环的出现时间建立了先后顺序的关系，容易想到把环缩成点，我们最终可以得到一张 DAG。
我们不难发现，只要连出来的有向图是 DAG（即没有环），就一定合法。
一张 DAG 对应若干方案，我们考虑对每一种 DAG 计算他的贡献，问题变为 DAG 计数，以及计算一张 DAG 的贡献。
DAG 计数是经典的，容斥即可。
对于一张 DAG，我们考虑每个环的后继环（包括自己）的总点数 \(p\)，枚举这个环的任意一点 \(i\) 的 \(a_i\) 在 \(b_i\) 中的位置 \(j\)，得到方案数为 \(A_{n-p}^{j-1}\times (n-j)!\)。
这个东西结合 DAG 计数，是容易的。
但是 \(n\le 40\)，我们不能直接状压。注意到我们只关心状压集合 \(S\) 中所有环的大小，把 \(S\) 中环的大小的可重集 \(S'\) 作为状压集合即可，事实证明状态数不大。
总结：
-
抽象模型，需要转化，常见转图论模型。
-
尝试从不同的角度出发，不要死磕
-
不要畏难，世界上没有不可做的题，只是不愿意思考罢了。
点击查看代码
#include<bits/stdc++.h>
#define ll long long
#define pir pair<ll,ll>
#define mkp make_pair
#define fi first
#define se second
#define pb push_back
using namespace std;
const ll maxn=42, mod=1e9+7;
ll n,a[maxn],cnt[maxn],d[maxn],siz[maxn],r[maxn],fac[maxn],ifac[maxn],f[10010];
ll find(ll x) {return x==d[x]? x:d[x]=find(d[x]);}
ll to[maxn];
ll Encode(){
ll S=0;
for(ll i=1;i<=n;i++) S=S*(cnt[i]+1)+to[i];
return S+1;
} ll tot=0;
ll power(ll a,ll b=mod-2){
ll s=1;
while(b){
if(b&1) s=s*a%mod;
a=a*a%mod; b>>=1;
} return s;
}
void dfs2(ll x,ll &ret){
if(x>n){
ll S=Encode();
if(S==tot) return;
ll c=mod-1, prod=1, p=0;
for(ll i=1;i<=n;i++) p+=r[i]*i,
prod=prod*fac[r[i]]%mod*ifac[to[i]]%mod*ifac[r[i]-to[i]]%mod;
for(ll i=1;i<=n;i++){
if(to[i]==r[i]) continue;
if((r[i]-to[i])&1) c=mod-c;
ll sum=0;
for(ll j=1;j<=n-p+1;j++)
sum=(sum+fac[n-p]*ifac[n-p-(j-1)]%mod*fac[n-j])%mod;
prod=prod*power(sum,(r[i]-to[i])*i)%mod;
} prod=prod*f[S]%mod;
ret=(ret+c*prod)%mod; return;
}
for(ll i=0;i<=r[x];i++)
to[x]=i, dfs2(x+1,ret);
}
void dfs(ll i){
if(i>n){ ++tot;
if(tot==1){
f[1]=1; return;
}
for(ll i=1;i<=n;i++) r[i]=to[i];
dfs2(1,f[tot]);
return;
}
for(ll j=0;j<=cnt[i];j++){
to[i]=j; dfs(i+1);
}
}
int main(){
scanf("%lld",&n); fac[0]=1;
for(ll i=1;i<=n;i++) fac[i]=fac[i-1]*i%mod;
ifac[n]=power(fac[n]);
for(ll i=n;i;i--) ifac[i-1]=ifac[i]*i%mod;
for(ll i=1;i<=n;i++) d[i]=i, siz[i]=1;
for(ll i=1;i<=n;i++){
scanf("%lld",a+i);
if(find(i)!=find(a[i])){
siz[find(a[i])]+=siz[find(i)];
d[find(i)]=find(a[i]);
}
}
for(ll i=1;i<=n;i++)
if(d[i]==i) ++cnt[siz[i]];
dfs(1);
printf("%lld",f[tot]);
return 0;
}