题意

输入一个长度为 $n$ 的数组 $g_i$，满足 $g_i \in [1,n]$。现在问有多少个长度为 $n$ 的数组 $f_i$，满足：

• $f_i \in [1,n]$
• $g_{f_i} = f_{g_i}$

$n \leq 5000$

题解

看到 $g$ 的定义，立马想到建图连边 $i \to g_i$，这样连出的是一个基环内向树森林。

可以先从特殊情况开始考虑。基环树的特殊情况无非就是树（某个点带有自环）和环。这题先从环开始入手比较方便。

假设我们现在知道了 $f_i,g_i$，由于 $f_{g_i} = g_{f_i}$，所以我们就可以推出 $g_{f_i}$ 是什么。所以枚举环上随便一个点的 $g$，然后推一遍就好了（注意最后要检查推回来的是否合法）

那么对于树，难点在于我们不能去枚举所有的叶子，所以我们考虑从上往下推。假设知道了 $f_v,g_v$，对于 $v$ 的一个孩子 $x$（也就是满足 $f_x=v$ 的点），有 $g_{f_x} = f_{g_x}$，所以就可以得到 $g_x$ 就是满足 $f_{g_x} = g_{f_x}$ 的所有点，也就是 $g_{f_x}$ 的儿子。所以设 $f_{v,i}$ 表示钦定 $g_v=i$ ，$v$ 子树的方案数。转移枚举每个儿子填哪个数就好了。复杂度均摊一下是 $O(n^2)$ 的。

对于基环树森林，首先每个弱联通块互相独立，先把每个小树 dp 完之后，再枚举环上的一个点，跑一遍加起来就好了。

分离基环树的方法：拓扑排序，标记每个点是否在环上，然后就可以找所有环上的非环邻点去 dfs 了。

代码

#include <bits/stdc++.h>

#define fi first
#define se second
#define DB double
#define U unsigned
#define P std::pair
#define LL long long
#define LD long double
#define pb emplace_back
#define MP std::make_pair
#define SZ(x) ((int)x.size())
#define all(x) x.begin(),x.end()
#define CLR(i,a) memset(i,a,sizeof(i))
#define FOR(i,a,b) for(int i = a;i <= b;++i)
#define ROF(i,a,b) for(int i = a;i >= b;--i)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl

const int MAXN = 5000+5;
const int ha = 1e9 + 7;
int to[MAXN],n,deg[MAXN];
int f[MAXN][MAXN];
std::vector<int> G[MAXN];
bool vis[MAXN];

inline void add(int &x,int y){
    x += y-ha;x += x>>31&ha;
}

int tmp[MAXN];

inline void dfs(int v){
    FOR(i,1,n) f[v][i] = 1;
    for(auto x:G[v]){
        if(deg[x]) continue;
        dfs(x);
        FOR(i,1,n){
            int s = 0;
            for(auto y:G[i]) add(s,f[x][y]);
            f[v][i] = 1ll*f[v][i]*s%ha;
        }
    }
}

int g[MAXN];

class CommutativeMapping {
public:
    inline int count(std::vector<int> vec){
        n = SZ(vec);FOR(i,1,n) to[i] = vec[i-1]+1,++deg[to[i]],G[to[i]].pb(i);
        std::queue<int> q;FOR(i,1,n) if(!deg[i]) q.push(i);
        while(!q.empty()){
            int v = q.front();q.pop();
            if(!--deg[to[v]]) q.push(to[v]);
        }
        FOR(i,1,n) if(deg[i]) dfs(i);
        int ans = 1;
        FOR(i,1,n){
            if(!deg[i] || vis[i]) continue;
            int v = i;
            std::vector<int> tmp;
            do{
                vis[v] = 1;tmp.pb(v);
                v = to[v];
            }while(!vis[v]);
            int sm = 0;
            FOR(j,1,n){
                g[0] = j;
                FOR(k,0,SZ(tmp)-1){
                    g[k+1] = to[g[k]];
                }
                if(g[SZ(tmp)] != g[0]) continue;
                int coe = 1;
                FOR(k,0,SZ(tmp)-1) coe = 1ll*coe*f[tmp[k]][g[k]]%ha;
                add(sm,coe);
            }
            ans = 1ll*ans*sm%ha;
        }
        return ans;
    }
};