考虑容斥，计算钦定 \(k\) 位满足 \(r_i=p_i\) 或 \(r_i=q_i\) 的方案数。
建出 \(n\) 个点，将每对 \(p_i,q_i\) 连边，由于每个点度数都是 \(2\)，所以会形成若干个环和一些孤立点(自环)。
对于自环，则有 \(p_i=q_i\)，只用考虑 \(r_i\) 是否等于 \(p_i\)，这可以用 DP 轻松解决。
不考虑自环，那么选出 \(k\) 个位置产生冲突就等价于在图上选择 \(k\) 条边，并且给每条被选择的边分配一个它的端点，并且每个点只能分配给一条边。
显然可以把图中的每个环分开考虑。
如果把环上的边也看成点，对于原先的一条边 \((u,v)\)，断开 \((u,v)\) 并新建点 \(w\)，连接 \((u,w),(w,v)\)，容易看出修改后的图仍然是环，此时选出 \(k\) 个位置就等价于选择 \(k\) 对匹配点，匹配的意思是选择相邻两个点。
考虑对于一个 \(n\) 个点的环，选择 \(k\) 对匹配的方案数是多少。
如果 \((1,n)\) 不匹配，那么就是先把 \(k\) 个点拿出，然后在剩下 \(n-k\) 个点中选出 \(k\) 个点作为匹配的右端点，最后把先前拿出的 \(k\) 个点插到选择的 \(k\) 个点的左边，这样就形成了 \(k\) 对匹配，则方案数为 \(\dbinom{n-k}{k}\)。
否则 \((1,n)\) 匹配的话，就是要在 \(2,\cdots,n-1\) 中选出 \(k-1\) 对匹配，则方案数为 \(\dbinom{n-2-(k-1)}{k-1}=\dbinom{n-k-1}{k-1}\)。
于是对于原问题中一个大小为 \(x\) 的环，选出 \(k\) 对匹配的方案数为 \(\dbinom{2\times x-k}{k}+\dbinom{2\times x-k-1}{k-1}\)。
用背包 \(\mathcal O(n^2)\) 合并每个环的方案数，计算出整张图钦定 \(k\) 位不合法的方案数 \(f_k\)，由于剩下的位置可以随便填，那么方案数是 \(f_k\times (n-k)!\)。
由容斥得最终答案是
\[\sum_{i=0}^n (-1)^i f_i\times (n-i)! \]Code：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 6005, mod = 1e9 + 7;
int n;
int a[N], b[N];
int fac[N], inv[N];
bool vis[N];
int f[N], g[N], tot;
int qpow(int x, int y) {
int res = 1;
while (y) {
if (y & 1) res = 1ll * res * x % mod;
x = 1ll * x * x % mod;
y >>= 1;
}
return res;
}
void init(int n) {
fac[0] = 1;
for (int i = 1; i <= n; ++i) fac[i] = 1ll * fac[i - 1] * i % mod;
inv[n] = qpow(fac[n], mod - 2);
for (int i = n - 1; ~i; --i) inv[i] = 1ll * inv[i + 1] * (i + 1) % mod;
}
int C(int n, int m) {
if (n < 0 || m < 0 || n < m) return 0;
return 1ll * fac[n] * inv[n - m] % mod * inv[m] % mod;
}
void add(int &a, int b) {
a += b;
if (a >= mod) a -= mod;
if (a < 0) a += mod;
}
int main() {
scanf("%d", &n);
for (int i = 1; i <= n; ++i) scanf("%d", &a[i]);
for (int i = 1; i <= n; ++i) scanf("%d", &b[a[i]]);
init(n * 2);
f[0] = 1;
for (int k = 1; k <= n; ++k) if (!vis[k]) {
int cnt = 0;
for (int i = k; !vis[i]; i = b[i]) vis[i] = 1, ++cnt;
memset(g, 0, sizeof g);
if (cnt == 1) {
g[0] = f[0];
for (int i = 1; i <= tot + 1; ++i) g[i] = (f[i - 1] + f[i]) % mod;
}
else {
for (int i = 0; i <= tot; ++i)
for (int j = 0; j <= cnt; ++j)
add(g[i + j], 1ll * f[i] * ((C(2 * cnt - j, j) + C(2 * cnt - j - 1, j - 1)) % mod) % mod);
}
memcpy(f, g, sizeof f); tot += cnt;
}
int ans = 0;
for (int i = 0; i <= n; ++i) {
f[i] = 1ll * f[i] * fac[n - i] % mod;
if (i & 1) add(ans, -f[i]);
else add(ans, f[i]);
}
printf("%d", ans);
return 0;
}