[HEOI2016/TJOI2016]求和 - 斯特林数 + NTT

Description

计算函数的值

\[f(n) = \sum \limits_{i=0}^{n} \sum \limits_{j=0}^{i} 2^j \times j! \times S(i,j) \]

Solution

大家好，我是练习推柿子半天的个人练习生\(newbielyx\)。

\[f(n) = \sum \limits_{i=0}^{n} \sum \limits_{j=0}^{i} 2^j \times j! \times S(i,j) \]

\[= \sum \limits_{i=0}^{n} \sum \limits_{j=0}^{i} 2^j \times j! \times \frac{1}{j!} \times \left (\sum \limits_{k=0}^{j} (-1)^k C(j,k)(j-k)^i\right) \]

\[= \sum \limits_{i=0}^{n} \sum \limits_{j=0}^{i} 2^j \sum \limits_{k=0}^{j} (-1)^k C(j,k)(j-k)^i \]

\[= \sum \limits_{j=0}^{n}2^j \sum \limits_{i=j}^{n} \sum \limits_{k=0}^{j} (-1)^k C(j,k) (j-k)^i \]

\[= \sum \limits_{j=0}^{n} 2^j \sum \limits_{k=0}^{j} (-1)^k C(j,k) \sum \limits_{i=0}^{n}(j-k)^i \]

\[= \sum \limits_{j=0}^{n} 2^j \sum \limits_{k=0}^{j} (-1)^k \frac{j!}{k!(j-k)!} \sum \limits_{i=0}^{n}(j-k)^i \]

\[= \sum \limits_{j=0}^{n} 2^j \cdot j! \sum \limits_{k=0}^{j} \frac{(-1)^k}{k!} \frac{\sum \limits_{i=0}^{n}(j-k)^i}{(j-k)!} \]

容易看出上面是一个等比数列，继续推倒

\[f(n) = \sum \limits_{j=0}^{n} 2^j \cdot j! \sum \limits_{k=0}^{j} \frac{(-1)^k}{k!} \frac{(j-k)^{n+1}-1}{(j-k-1)(j-k)!} \]

\[= \sum \limits_{j=0}^{n} 2^j \cdot j! \sum \limits_{k=0}^{j} \frac{(-1)^k}{k!} \frac{(j-k)^{n+1}-1}{(j-k)!(j-k-1)} \]

然后你就可以发现，卷起来了！

Code

#include <bits/stdc++.h>
using namespace std;

inline int ty() {
	char ch = getchar(); int x = 0, f = 1;
	while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
	while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
	return x * f;
}

typedef long long ll;
const ll P = 998244353, G = 3, Gx = 332748118;
const int _ = 4e5 + 10;
int N, r[_];
ll fac[_], facinv[_], A[_], B[_];

ll fpow(ll a, ll b) {
	ll ret = 1;
	for ( ; b; b >>= 1) {
		if (b & 1) ret = ret * a % P;
		a = a * a % P;
	}
	return ret;
}

inline void pre() {
	fac[0] = fac[1] = 1;
	for (int i = 2; i <= N; ++i) fac[i] = fac[i - 1] * i % P;
	facinv[0] = 1, facinv[N] = fpow(fac[N], P - 2);
	for (int i = N - 1; i >= 1; --i) facinv[i] = facinv[i + 1] * (i + 1) % P;
	for (int i = 0; i <= N; ++i) {
		A[i] = (i & 1) ? -1 : 1;
		A[i] = (A[i] * facinv[i] + P) % P;
	}
	for (int i = 2; i <= N; ++i) {
		B[i] = (fpow(i, N + 1) - 1 + P) % P;
		B[i] = B[i] * facinv[i] % P;
		B[i] = B[i] * fpow(i - 1, P - 2) % P;
	}
	B[0] = 1, B[1] = N + 1;
}

void NTT(int lim, ll *a, int op) {
	for (int i = 0; i < lim; ++i)
		if (i < r[i]) swap(a[i], a[r[i]]);
	for (int len = 2; len <= lim; len <<= 1) {
		int mid = len >> 1;
		ll Wn = fpow(op == 1 ? G : Gx, (P - 1) / len);
		for (int i = 0; i < lim; i += len) {
			ll w = 1;
			for (int j = 0; j < mid; ++j, w = w * Wn % P) {
				ll x = a[i + j], y = w * a[i + j + mid] % P;
				a[i + j] = (x + y) % P;
				a[i + j + mid] = (x - y + P) % P;
			}
		}
	}
}

inline void solve() {
	int lim = 1, k = 0;
	while (lim <= N + N) lim <<= 1, ++k;
	for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (k - 1));
	NTT(lim, A, 1);
	NTT(lim, B, 1);
	for (int i = 0; i < lim; ++i) A[i] = (A[i] * B[i]) % P;
	NTT(lim, A, -1);
	ll inv = fpow(lim, P - 2);
	for (int i = 0; i <= N + N; ++i) A[i] = A[i] * inv % P;
}

int main() {
#ifndef ONLINE_JUDGE
	freopen("sum.in", "r", stdin);
	freopen("sum.out", "w", stdout);
#endif
	N = ty();
	pre();
	solve();
	ll ans = 0, t = 1;
	for (int i = 0; i <= N; ++i, t = t * 2 % P) ans = (ans + t * fac[i] % P * A[i] % P) % P;
	printf("%lld\n", ans);
	return 0;
}