BZOJ4555[HEOI2016/TJOI2016]求和

题目链接

洛谷

BZOJ

前置知识

第二类斯特林数

含义

(S_n^m)表示将(n)个互不相同的元素划分为(m)个非空集合的方案数

递推式

[S_n^m = S_{n - 1}^{m - 1} + m cdot S_{n - 1}^m ]

通项公式

[S_n^m = frac{1}{m!} cdot sum_{k = 0}^m (-1)^k cdot {m choose k} cdot (m - k)^n ]

NTT

类似于(FFT),在固定模数的情况下,可以将单位根用模数(P)的原根的(frac{P - 1}{n})次幂代替,避免复数运算掉精度

还是放个链接吧:NTT详解

解析

首先(j)的上界可以扩展到(n),因为(j > i)的时候(S(i,j) = 0)

后两项都只和(j)有关,所以换一下枚举顺序,把这两项提出来:

[f(n) = sum_{j = 0}^n 2^j cdot (j!) cdot sum_{i = 0}^n S(i, j) ]

你看那个(S(i,j))就很不友好,我们直接把它干掉,用通项替换:

[f(n) = sum_{j = 0}^n 2^j cdot (j!) cdot sum_{i = 0}^n frac{1}{j!} cdot sum_{k = 0}^j (-1)^k cdot frac{j!}{k! cdot (j - k)!} cdot (j - k)^i ]

然后一波化简可得:

[f(n) = sum_{j = 0}^n 2^j cdot (j!) cdot sum_{i = 0}^n sum_{k = 0}^j (-1)^k cdot frac{1}{k! cdot (j - k)!} cdot (j - k)^i ]

我们把后面那一坨提出来,令:

[g(j) = sum_{i = 0}^n sum_{k = 0}^j (-1)^k cdot frac{1}{k! cdot (j - k)!} cdot (j - k)^i \ = sum_{k = 0}^j frac{(-1)^k}{k!} cdot frac{sum_{i = 0}^n (j - k)^i}{(j - k)!} ]

再令:

[a(i) = frac{(-1)^i}{i!},b(i) = frac{sum_{k = 0}^n i^k}{i!} = frac{i^{k + 1} - 1}{(i - 1) cdot i!} ]

容易发现(g)就是(a)(b)的卷积

特别地,(b(0) = 1, b(1) = n + 1)

那么用(NTT)或者(FFT)什么的求出(g)后,答案就是(f(n) = sum_{i = 0}^n 2^i cdot (i!) cdot g(i))

代码

我写的(NTT),然而忘了开(4n)的数组(RE)(3)发。。。。。

#include <cstdio>
#include <cstring>
#include <iostream>
#define MAXN 100005

typedef long long LL;
const LL mod = 998244353ll;
LL N, g[MAXN << 2], a[MAXN << 2], b[MAXN << 2], fact[MAXN];
LL G = 3, ans;

LL qpower(LL, LL);
void NTT(LL *, int, int);
inline LL inverse(LL x) { return qpower(x, mod - 2); }
int main() {
	std::ios::sync_with_stdio(false);
	std::cin >> N;
	fact[0] = 1;
	for (int i = 1; i <= N; ++i) fact[i] = fact[i - 1] * i % mod;
	for (int i = 0; i <= N; ++i) a[i] = ((i & 1) ? -1 : 1) * inverse(fact[i]);
	b[0] = 1, b[1] = N + 1;
	for (int i = 2; i <= N; ++i) b[i] = (qpower(i, N + 1) - 1) * inverse(fact[i]) % mod * inverse(i - 1) % mod;
	int sz = 1;
	while (sz <= (N << 1)) sz <<= 1;
	NTT(a, sz, 1);
	NTT(b, sz, 1);
	for (int i = 0; i < sz; ++i) g[i] = a[i] * b[i] % mod;
	NTT(g, sz, -1);
	LL pow2 = 1;
	for (int i = 0; i <= N; ++i, pow2 = (pow2 << 1) % mod)
		ans = (ans + pow2 * fact[i] % mod * g[i] % mod) % mod;
	std::cout << ans << std::endl;
	return 0;
}
LL qpower(LL x, LL y) {
	LL res = 1;
	while (y) {
		if (y & 1) (res *= x) %= mod;
		(x *= x) %= mod;
		y >>= 1;
	}
	return res;
}
void NTT(LL *arr, int sz, int type) {
	int rev[MAXN << 2]; rev[0] = 0;
	for (int i = 1; i < sz; ++i) {
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (sz >> 1) : 0);
		if (rev[i] > i) std::swap(arr[i], arr[rev[i]]);
	}
	for (int len = 2, half = 1; len <= sz; len <<= 1, half <<= 1) {
		LL wn = qpower(G, (mod - 1) / len);
		if (type == -1) wn = inverse(wn);
		for (int i = 0; i < sz; i += len) {
			LL w = 1;
			for (int j = 0; j < half; ++j, w = w * wn % mod) {
				LL tmp1 = arr[i + j], tmp2 = arr[i + half + j] * w % mod;
				arr[i + j] = tmp1 + tmp2; if (arr[i + j] >= mod) arr[i + j] -= mod;
				arr[i + half + j] = tmp1 - tmp2; if (arr[i + half + j] < 0) arr[i + half + j] += mod;
			}
		}
	}
	if (type == -1) {
		LL inv = inverse(sz);
		for (int i = 0; i < sz; ++i) arr[i] = arr[i] * inv % mod;
	}
}
//Rhein_E
原文地址:https://www.cnblogs.com/Rhein-E/p/10439527.html