题意
给定长度为 $N$ 的两个排列 $p$ 和 $q$，求长度为 $N$ 的排列 $r$ 的个数，使得 $r_i \neq p_i$ 且 $r_i \neq q_i$，对于任意的 $1 \leq i \leq N$ 都成立。
$1 \leq N \leq 3000$
题解
设 $f_i$ 表示任意选 $i$ 个位置，对于这些位置都有 $r_i = p_i$ 或 $r_i = q_i$。根据容斥，答案是
$$\sum\limits_{0 \leq i \leq N}\left(-1\right)^if_i\left(n - i\right)!$$
考虑 $f_i$ 的计算。由于 $r_i$ 不能重复导致不好计数，我们可以换一种思路，让限制匹配位置。具体来说，在 $p_i$ 和 $q_i$ 之间连边，如此在图中，限制 $i$ 被满足表示为把 $p_i$ 或 $q_i$ 分配给了 $\left(p_i, q_i\right)$ 这条边。
还是不好计数，但还有救，把 $\left(p_i, q_i\right)$ 变成 $\left(p_i, s_i\right)$ 和 $\left(s_i, q_i\right)$，相当于在这两条边之间新增一个点。于是点到边的分配被转化成了点到点的匹配。
没什么用，我们需要更多性质。发现图中每个点刚好有一个入和一个出，构成若干个环（可能包含自环）。环之间显然是独立的。于是我们可以对环求出答案然后 DP 合并。
环内的答案是好计算的。假设有 $n$ 个点（包含我们新增的 $s$ 系列点），我们要选出 $k$ 组匹配。如果 $\left(1, n\right)$ 不被选为匹配，那么方案数是 $\binom{n - k}k$。如果 $\left(1, n\right)$ 被选为匹配，那么方案数是 $\binom{n - 2 - \left(k - 1\right)}{k - 1}$。于是相加得到总方案数
$$\binom{n - k}k + \binom{n - k - 1}{k - 1}$$
其实是有 $\tilde{\mathcal O}\left(n\right)$ 做法，不过要多项式，果断放弃。想学可以看这里。
代码
/* Please submit with C++17! It's best to use C++20 or higher version.
* No header file and no RBLIB (https://git.rbtr.ee/root/Template).
* By Koicy (https://koicy.ly)
* Email n@rbtr.ee
* I've reached the end of my fantasy.
__ __ __ _
_____/ /_ / /_________ ___ / /______ (_)______ __
/ ___/ __ \/ __/ ___/ _ \/ _ \ ______ / //_/ __ \/ / ___/ / / /
/ / / /_/ / /_/ / / __/ __/ /_____/ / ,< / /_/ / / /__/ /_/ /
/_/ /_.___/\__/_/ \___/\___/ /_/|_|\____/_/\___/\__, /
SIGN /___*/
#ifndef XCODE
constexpr bool _CONSOLE = false;
#else
constexpr bool _CONSOLE = true;
#endif
#define __NO_MAIN__ false
#define __ENABLE_RBLIB__ true
constexpr bool _MTS = false, SPC_MTS = false;
constexpr char EFILE[] = "";
#define FULL(arg) arg.begin(), arg.end()
// :/
signed STRUGGLING([[maybe_unused]] unsigned long TEST_NUMBER) {
constexpr tp mod = 1e9 + 7;
tp n = bin;
vetp p = bin.vcr(n, 1);
vetp q = bin.vcr(n, 1);
vector<vetp> e(n + 1);
for (tp i = 1; i <= n; ++i) {
e[p[i]].push_back(q[i]);
e[q[i]].push_back(p[i]);
}
vector<tp> sec { 0 };
set<tp> vis;
for (tp i = 1; i <= n; ++i) {
auto dfs = [&](auto $, tp x, tp& c) -> void {
if (vis.count(x)) return;
++c;
vis.insert(x);
for (auto& j : e[x]) $($, j, c);
};
if (vis.count(i)) continue;
sec.emplace_back();
dfs(dfs, i, sec.back());
}
vector<vetp> f(sec.size(), vetp(n + 1));
vetp fact(2 * n + 1, 1);
vetp ifac(2 * n + 1);
for (tp i = 2; i <= 2 * n; ++i) fact[i] = fact[i - 1] * i % mod;
for (tp i = 0; i <= 2 * n; ++i) ifac[i] = RCAL::inverse(fact[i], mod);
auto binom = [&](tp n, tp m) -> tp {
if (n < 0 || m < 0 || n < m) return 0;
return fact[n] * ifac[m] % mod * ifac[n - m] % mod;
};
f[0][0] = 1;
for (tp i = 1; i < sec.size(); ++i) {
if (sec[i] == 1) {
f[i][0] = f[i - 1][0];
for (tp j = 1; j <= n; ++j) f[i][j] = (f[i - 1][j] + f[i - 1][j - 1]) % mod;
continue;
}
for (tp j = 0; j <= n; ++j) {
for (tp k = 0; k <= sec[i] && j + k <= n; ++k) {
tp l = 2 * sec[i];
f[i][j + k] = (f[i][j + k] + f[i - 1][j] * (binom(l - k, k) + binom(l - k - 1, k - 1))) % mod;
}
}
}
tp tar = 0;
for (tp i = 0; i <= n; ++i) {
tp t = (i & 1 ? -1 : 1) * fact[n - i];
tar = (tar + t * f.back()[i]) % mod;
}
bin << (tar + mod) % mod << '\n';
return 0;
}
void MIST() {
}
// :\ */