集合划分个数

题目链接

集合划分个数

做法

就是求 $ B_n $（第 $ n $ 项贝尔数，将 $ n $ 个有标号的球划分为若干集合的方案数）。
一个非空集合的指数生成函数为 $ F(x) = e^x - 1 $ ，则：

\[B(x) = \sum_{i = 0}^{\infin} \frac{F^i(x)}{i!}\\\\ B(x) = \exp(F(x)) = \exp(e^x-1) \]

#include <bits/stdc++.h>
#define pb push_back
using namespace std;
typedef long long ll;
const int mod = 998244353;

inline int Max(const int &x, const int &y) { return x > y ? x : y; }
inline int Min(const int &x, const int &y) { return x < y ? x : y; }
inline int add(const int &x, const int &y) { return x + y < mod ? x + y : x + y - mod; }
inline int sub(const int &x, const int &y) { return x - y < 0 ? x - y + mod : x - y; }
inline int mul(const int &x, const int &y) { return (int)((ll)x * y % mod); }
inline int ksm(int x, int y = mod - 2) {
	int ss = 1; for(; y; y >>= 1, x = mul(x, x)) if(y & 1) ss = mul(ss, x); return ss;
}
namespace Poly {
	inline int Get(int x) { int ss = 1; for(; ss <= x; ss <<= 1); return ss; }
	void ntt(vector<int> &A, int lmt, int opt) {
		A.resize(lmt + 5);
		for(int i = 0, j = 0; i < lmt; i++) {
			if(i > j) swap(A[i], A[j]);
			for(int k = lmt >> 1; (j ^= k) < k; k >>= 1);
		}
		vector<int> w(lmt >> 1);
		for(int mid = 1; mid < lmt; mid <<= 1) {
			w[0] = 1;
			int w0 = ksm(opt == 1 ? 3 : (mod + 1) / 3, (mod - 1) / mid / 2);
			for(int j = 1; j < mid; j++) w[j] = mul(w[j - 1], w0);
			for(int R = mid << 1, j = 0; j < lmt; j += R)
				for(int k = 0; k < mid; k++) {
					int x = A[j + k], y = mul(w[k], A[j + mid + k]);
					A[j + k] = add(x, y), A[j + mid + k] = sub(x, y);
				}
		}
		if(opt == -1)
			for(int i = 0, inv = ksm(lmt); i < lmt; i++) A[i] = mul(A[i], inv);
	}
	vector<int> Add(const vector<int> &a, const vector<int> &b) {
		vector<int> res(Max(a.size(), b.size()));
		for(int i = 0; i < (int)res.size(); i++) {
			if(i < (int)a.size()) res[i] = add(res[i], a[i]);
			if(i < (int)b.size()) res[i] = add(res[i], b[i]);
		}
		return res;
	}
	vector<int> Mul(const vector<int> &a, const vector<int> &b) {
		vector<int> A = a, B = b; int lmt = Get(a.size() + b.size() - 2);
		ntt(A, lmt, 1), ntt(B, lmt, 1);
		for(int i = 0; i < lmt; i++) A[i] = mul(A[i], B[i]);
		ntt(A, lmt, -1); return A.resize(a.size() + b.size() - 1), A;
	}
	vector<int> Inv(const vector<int> &A, int sz = -1) {
		if(sz == -1) sz = A.size();
		vector<int> res; if(sz == 1) return res.pb(ksm(A[0])), res;
		res = Inv(A, (sz + 1) / 2);
		vector<int> tmp(A.begin(), A.begin() + sz);
		int lmt = Get(sz * 2 - 2);
		ntt(tmp, lmt, 1), ntt(res, lmt, 1);
		for(int i = 0; i < lmt; i++) res[i] = mul(sub(2, mul(res[i], tmp[i])), res[i]);
		ntt(res, lmt, -1); return res.resize(sz), res;
	}
	void Div(const vector<int> &A, const vector<int> &B, vector<int> &D, vector<int> &R) {
		if(B.size() > A.size()) return (void)(D.clear(), D.pb(0), R = A);
		vector<int> a = A, b = B, iB; int n = A.size(), m = B.size();
		reverse(a.begin(), a.end()), reverse(b.begin(), b.end());
		b.resize(n - m + 1), iB = Inv(b, n - m + 1);
		D = Mul(a, iB), D.resize(n - m + 1), reverse(D.begin(), D.end());
		R = Mul(B, D);
		for(int i = 0; i < m - 1; i++) R[i] = (mod + A[i] - R[i]) % mod;
		R.resize(m - 1);
	}
	vector<int> Ln(const vector<int> &A) {
		vector<int> f, g = Inv(A); f.resize(A.size());
		for(int i = 1; i < (int)f.size(); i++) f[i - 1] = mul(i, A[i]);
		f[f.size() - 1] = 0, f = Mul(f, g), f.resize(A.size());
		for(int i = (int)f.size() - 1; i > 0; i--) f[i] = mul(f[i - 1], ksm(i));
		f[0] = 0; return f;
	}
	vector<int> Exp(const vector<int> &A, int sz = -1) {
		if(sz == -1) sz = A.size();
		vector<int> res; if(sz == 1) return res.pb(1), res;
		res = Exp(A, (sz + 1) / 2), res.resize(sz); vector<int> tmp = Ln(res);
		for(int i = 0; i < (int)tmp.size(); i++) tmp[i] = add(sub(0, tmp[i]), A[i]);
		tmp[0] = add(1, tmp[0]), res = Mul(res, tmp);
		return res.resize(sz), res;
	}
	vector<int> Ksm(const vector<int> &A, const int &K) {
		int sz = A.size(), low = 0; vector<int> ta, res;
		for(int i = 0; i < sz; i++) if(A[i]) { low = i; break; }
		for(int i = low; i < sz; i++) ta.pb(A[i]);
		int mu = ta[0], inv = ksm(mu), Mu = ksm(mu, K);
		for(int i = 0; i < (int)ta.size(); i++) ta[i] = 1ll * ta[i] * inv % mod;
		ta = Ln(ta);
		for(int i = 0; i < (int)ta.size(); i++) ta[i] = 1ll * ta[i] * K % mod;
		ta = Exp(ta);
		for(int i = 0; i < (int)ta.size(); i++) ta[i] = 1ll * ta[i] * Mu % mod;
		int cnt = 0;
		for(int i = 0; i < low; i++)
			for(int j = 1; j <= K; j++) { ++cnt, res.pb(0); if(cnt > sz) break; }
		for(int i = 0; i < (int)ta.size(); i++) res.pb(ta[i]);
		return res.resize(sz), res;
	}
	vector<int> Ksm(const vector<int> &A, const vector<int> &K) {
		int mu = 0;
		for(int i = 0; i < (int)K.size(); i++) mu = add(mul(mu, 10), K[i]);
		return Ksm(A, mu);
	}
	vector<int> Sqrt(const vector<int> &A, int sz = -1) {
		if(sz == -1) sz = A.size();
		vector<int> res; if(sz == 1) { return res.pb(1), res; }
		res = Sqrt(A, (sz + 1) / 2); vector<int> tmp(A.begin(), A.begin() + sz);
		res.resize(sz), tmp = Mul(tmp, Inv(res));
		for(int i = 0; i < sz; i++) res[i] = 1ll * (res[i] + tmp[i]) * (mod + 1) / 2 % mod;
		return res;
	}
	vector<int> tt[4000010], res;
	vector<int> Evaluate_init(int u, int l, int r, const vector<int> &a) {
		vector<int> ttt;
		if(l >= r) { ttt.pb(mod - a[l]), ttt.pb(1); return tt[u] = ttt; }
		int mid = (l + r) >> 1;
		Evaluate_init(u * 2, l, mid, a), Evaluate_init(u * 2 + 1, mid + 1, r, a);
		return tt[u] = Mul(tt[u * 2], tt[u * 2 + 1]);
	}
	void Evaluate(int u, int l, int r, const vector<int> &f, const vector<int> &a) {
		if(r - l + 1 <= 512) {
			for(int k = l; k <= r; k++) {
				int w = 1, x = 0;
				for(int i = 0; i < (int)f.size(); i++) x = add(x, mul(w, f[i])), w = mul(w, a[k]);
				res.pb(x);
			}
			return ;
		}
		int mid = (l + r) >> 1; vector<int> tmp; Div(f, tt[u], tmp, tmp);
		Evaluate(u * 2, l, mid, tmp, a), Evaluate(u * 2 + 1, mid + 1, r, tmp, a);
	}
	vector<int> Evaluation(const vector<int> &f, const vector<int> &a) {
		Evaluate_init(1, 0, a.size() - 1, a);
		res.clear(), Evaluate(1, 0, a.size() - 1, f, a); return res;
	}
	
}
const int N = 100010;
int T;
vector<int> B;
int fac[N];

int main() {
	fac[0] = 1;
	for(int i = 1; i < 100010; i++) fac[i] = 1ll * fac[i - 1] * i % mod;
	B.resize(100010);
	for(int i = 1; i < (int)B.size(); i++) B[i] = ksm(fac[i]);
	B = Poly::Exp(B);
	scanf("%d", &T);
	for(int n; T; --T) {
		scanf("%d", &n);
		printf("%lld\n", 1ll * B[n] * fac[n] % mod);
	}
	return 0;
}