题意：很简单了，不再赘述。
做法：
直接做很困难，考虑容斥，钦定若干个 \(i=x_i\) 或者 \(i=y_i\)，然后如果钦定了 \(k\) 个，那么贡献是 \((n-k)!(-1)^k\)。
但是有个问题，我们不能无脑钦定，这些钦定条件间会有一些问题，比如我可能钦定 \(i=k\) 但是同时 \(j=k\)，这样就爆炸了。
所以我们考虑对于不能同时钦定的限制去连一条边，手玩一下你会很惊奇地发现这样会形成若干个环，并且这个环上不能取相邻元素！
具体怎么连环呢，对于同一个 \(i\) 间的两个是不能同存的，然后两行同一个数不能被同时钦定。
那么我们就可以对不同长度的环进行 dp，计算出来这个环选 \(x\) 个的有多少种方案，直接枚举第一个选或不选然后跑 dp 即可，复杂度 \(O(n^2)\)。
然后最后我们再用一次 dp 对所有的环统计出来钦定了 \(x\) 个限制的方案数，虽然是三重循环但是枚举第 \(x\) 个环再枚举环上选了多少个这两个是 \(O(n)\) 的。
最后乘上贡献即可，复杂度 \(O(n^2)\)。
注意一个细节，对于 \(x_i = y_i\) 的情况限制可以同时存在，这种需要特判一下，环上选一个的方案数只有 \(1\) 种。
代码：
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int maxn = 6005, mod = 1e9 + 7;
int f[maxn][maxn], t[2][maxn][2], use[maxn], dp[maxn][maxn], len[maxn], tot, jc[maxn];
int x[2][maxn], n;
void solve(int x) {memset(t, 0, sizeof(t));t[0][0][0] = 1;int cur = 0;for (int i = 2; i <= x; i++) {cur ^= 1;for (int j = 0; j <= (i + 1) / 2; j++)t[cur][j][0] = t[cur][j][1] = 0;for (int j = 0; j <= (i + 1) / 2; j++) {t[cur][j][0] = (t[cur ^ 1][j][0] + t[cur ^ 1][j][1]) % mod;t[cur][j][1] = t[cur ^ 1][j - 1][0];} }for (int i = 0; i <= x / 2; i++)f[x][i] = (t[cur][i][0] + t[cur][i][1]) % mod;memset(t, 0, sizeof(t));t[0][1][1] = 1;cur = 0;for (int i = 2; i <= x; i++) {cur ^= 1;for (int j = 0; j <= (i + 1) / 2; j++)t[cur][j][0] = t[cur][j][1] = 0;for (int j = 0; j <= (i + 1) / 2; j++) {t[cur][j][0] = (t[cur ^ 1][j][0] + t[cur ^ 1][j][1]) % mod;t[cur][j][1] = t[cur ^ 1][j - 1][0];} }for (int i = 0; i <= x / 2; i++)f[x][i] = (t[cur][i][0] + f[x][i]) % mod;
}
void solve() {for (int i = 1; i <= tot; i++)use[len[i]] = 1;for (int i = 1; i <= 2 * n; i++)if(use[i])solve(i);f[2][1] = 1;dp[0][0] = 1;int s = 0;for (int i = 1; i <= tot; i++) {s += len[i];for (int j = 0; j <= s / 2; j++)for (int k = 0; k <= min(len[i] / 2, j); k++)dp[i][j] = (dp[i][j] + dp[i - 1][j - k] * f[len[i]][k]) % mod;}int ans = 0;jc[0] = 1;for (int i = 1; i <= n; i++)jc[i] = jc[i - 1] * i % mod;for (int i = 0; i <= n; i++)ans = (ans + dp[tot][i] * jc[n - i] % mod * (i % 2 ? mod - 1 : 1) % mod) % mod;cout << ans << endl;
}
int pos[2][maxn], p[maxn], vis[maxn];
void dfs(int u, int d) {if(vis[u]) {len[tot] = d - vis[u];return ;}vis[u] = d;dfs(p[u], d + 1);
}
signed main() {cin >> n;for (int t = 0; t <= 1; t++) {for (int i = 1; i <= n; i++)cin >> x[t][i], pos[t][x[t][i]] = i;}for (int i = 1; i <= n; i++)p[i] = i + n, p[pos[1][x[0][i]] + n] = i;for (int i = 1; i <= 2 * n; i++)if(!vis[i]) tot++, dfs(i, 1);solve();return 0;
}
