BZOJ5406: Gift 第一类斯特林数

Description
定义两个排列相似度为一个排列交换两个元素得到另一个的最小步数。
给你两个排列$A,B$，其中一些元素是$0$，你可以补上一些数。
现在询问对于每一个$i$，补全后相似度为i的方案数。

---

Sample Input
3
1 0 0
0 2 0

---

Sample Output
1 2 1

---

我们可以知道对于两个排列，相似度其实就是$n-$环数。

根据两个序列可以连出一些边。

首先对于$x-x$这种边已经固定了，它要么已经形成了环，要么就可以缩起来。
易得一个$x-0$和$0-x$是无法直接构成一个环的，可以根据定义得知。

那么我们就将$x-0$和$0-x$分开考虑，最后再卷起来。
假设有$n$条$x-0$，有$k$条$0-0$，定义$f_i$为至少有$i$个环的方案数，可以得到下式：
$$f_i=\sum _{j=i}^n{C}_n^j{S}_j^i(n-j+k{)}^{\underline{n-j}}$$
$S$是第一类斯特林数
是这样理解的：就相当于你选$j$条$x-0$边强制构成$i$个环，然后剩下的$x-0$要么就是自由的组成环，要么就是将一些$0-0$边给缩起来。
然后就可以直接二项式反演得出有$i$个环的方案数。

直接把$x-0$和$0-x$卷起来之后，然后就剩下$0-0$的边，$0-0$的方案就相当于斯特林数乘个阶乘，直接卷起来就好了。

---

#include <ctime>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>

using namespace std;
typedef long long LL;
int _max(int x, int y) {return x > y ? x : y;}
int _min(int x, int y) {return x < y ? x : y;}
const LL mod = 998244353;
const int N = 2001;
int read() {
	int s = 0, f = 1; char ch = getchar();
	while(ch < '0' || ch > '9') {if(ch == '-') f = -1; ch = getchar();}
	while(ch >= '0' && ch <= '9') s = s * 10 + ch - '0', ch = getchar();
	return s * f;
}
void put(int x) {
	if(x >= 10) put(x / 10);
	putchar(x % 10 + '0');
}

int n, cnt1, cnt2, cnt3, cnt4;
int a[N], b[N], nxt[N * 2], ru[N * 2];
int C[N][N], S[N][N], A[N][N];
int f[N * 4], g[N * 4], h[N * 4], ans[N * 4];
int w1[N * 4], w2[N * 4];
bool v[N * 2];
int R[N * 4];

void add(int &x, int y) {
	y -= y >= mod ? mod : 0;
	x += y;
	x -= x >= mod ? mod : 0;
}

void dfs(int x, int tp) {
	v[x] = 1;
	if(nxt[x]) {
		if(v[nxt[x]]) cnt4++;
		else dfs(nxt[x], tp);
	} else {
		if(tp <= n && x > n) cnt1++;
		else if(tp > n && x <= n) cnt2++;
		else if(tp > n && x > n) cnt3++;
	}
}

int pow_mod(int a, int k) {
	int ans = 1;
	while(k) {
		if(k & 1) ans = (LL)ans * a % mod;
		a = (LL)a * a % mod; k /= 2;
	} return ans;
}

void pre(int len) {
	for(int i = 1; i < len; i <<= 1) {
		w1[i] = pow_mod(3, (mod - 1) / (i << 1));
		w2[i] = pow_mod(w1[i], mod - 2);
	}
}

void NTT(int y[], int len, int on) {
	for(int i = 0; i < len; i++) if(i < R[i]) swap(y[i], y[R[i]]);
	for(int i = 1; i < len; i <<= 1) {
		int wn = on == -1 ? w2[i] : w1[i];
		for(int j = 0; j < len; j += i << 1) {
			int w = 1;
			for(int k = 0; k < i; k++, w = (LL)w * wn % mod) {
				int u = y[j + k], v = (LL)y[j + k + i] * w % mod;
				y[j + k] = u, add(y[j + k], v);
				y[j + k + i] = u, add(y[j + k + i], mod - v);
			}
		} 
	} if(on == -1) {
		int inv = pow_mod(len, mod - 2);
		for(int i = 0; i < len; i++) y[i] = (LL)y[i] * inv % mod;
	}
}

void gao(int f[], int n) {
	for(int i = 0; i <= n; i++) {
		for(int j = i; j <= n; j++) {
			add(f[i], (LL)C[n][j] * S[j][i] % mod * A[n - j + cnt3][n - j] % mod);
		}
	} for(int i = 0; i <= n; i++) {
		int sum = 0;
		for(int j = i, d = 1; j <= n; j++, d = mod - d) {
			add(sum, (LL)C[j][i] * d % mod * f[j] % mod);
		} f[i] = sum;
	}
}

int main() {
	n = read();
	for(int i = 1; i <= n; i++) a[i] = read();
	for(int i = 1; i <= n; i++) b[i] = read();
	for(int i = 1; i <= 2 * n; i++) v[i] = 1;
	for(int i = 1; i <= n; i++) {
		if(!a[i]) a[i] = i + n; if(!b[i]) b[i] = i + n;
		v[a[i]] = v[b[i]] = 0;
		if(a[i] <= n || b[i] <= n) nxt[a[i]] = b[i], ++ru[b[i]];
	} for(int i = 1; i <= n * 2; i++) if(!v[i] && !ru[i]) dfs(i, i);
	for(int i = 1; i <= n * 2; i++) if(!v[i]) dfs(i, i);
	C[0][0] = A[0][0] = S[0][0] = 1;
	for(int i = 1; i <= n; i++) {
		A[i][0] = C[i][0] = 1;
		for(int j = 1; j <= i; j++) {
			C[i][j] = C[i - 1][j - 1], add(C[i][j], C[i - 1][j]);
			S[i][j] = S[i - 1][j - 1], add(S[i][j], (LL)S[i - 1][j] * (i - 1) % mod);
			A[i][j] = (LL)A[i][j - 1] * (i - j + 1) % mod;
		}
	} gao(f, cnt1), gao(g, cnt2);
	int len; for(len = 1; len <= n * 2; len <<= 1);
	for(int i = 0; i < len; i++) R[i] = R[i >> 1] >> 1 | (i & 1) * (len >> 1);
	pre(len), NTT(f, len, 1), NTT(g, len, 1);
	for(int i = 0; i < len; i++) h[i] = (LL)f[i] * g[i] % mod;
	for(int i = 0; i <= n; i++) f[i] = S[cnt3][i];
	for(int i = n + 1; i < len; i++) f[i] = 0;
	NTT(f, len, 1);
	for(int i = 0; i < len; i++) ans[i] = (LL)h[i] * f[i] % mod;
	NTT(ans, len, -1);
	for(int i = 0; i <= n; i++) ans[i] = (LL)ans[i] * A[cnt3][cnt3] % mod;
	for(int i = 0; i < n; i++) put(n - i - cnt4 >= 0 ? ans[n - i - cnt4] : 0), putchar(' ');
	return 0;
}