P4091 [HEOI2016/TJOI2016] 求和

P4091 [HEOI2016/TJOI2016] 求和

令 \(S(n,m)\) 表示第二类斯特林数。

对以下式子求和：

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

\(n\le 10^5\)（但是有线性做法

---

\(\forall m>n,S(n,m)=0\)，所以将第二个求和上界改为 \(n\)，答案不变。

考虑对于第二类斯特林数的容斥：

\[S(n,m)=\frac 1{m!}\sum_{k=0}^m(-1)^{m-k}\binom mkk^n \]

代入原式：

\[\begin{align*} &\sum_{i=0}^n\sum_{j=0}^i 2^j\times j!\times S(i,j) \\=&\sum_{i=0}^n\sum_{j=0}^n2^j\sum_{k=0}^j\binom jk(-1)^{j-k}k^i \\=&\sum_{j=0}^n2^j\sum_{k=0}^j\binom jk(-1)^{j-k}\sum_{i=0}^nk^i \end{align*} \]

令 \(a_k=\sum_{i=0}^nk^i\)，\(a_k\) 可以做到线性预处理。

\[\begin{align*} &\sum_{j=0}^n2^j\sum_{k=0}^j\binom jk(-1)^{j-k}a_k \\=&\sum_{j=0}^n\sum_{k=0}^j2^k2^{j-k}(-1)^{j-k}\frac {j!}{k!(j-k)!}a_k \\=&\sum_{0\le i+k\le n}(i+k)!\frac {(-2)^i}{i!}\frac{2^ka_k}{k!} \end{align*} \]

最后一步做了换元 \(i=j-k\)。可以发现这是一个卷积的形式，于是可以做到 \(\mathcal O(n\log n)\)。

---

能不能做到更好？

假如我们有一个生成函数 \(F(x)=\sum_{i=0}^na_ix^i\)。

那么 \(F(x+p)\) 中第 \(j\) 项的系数是多少？

\[\begin{align*} &[x^k]F(x+p) \\=&[x^k]\sum_{i=0}^na_i(x+p)^i \\=&[x^k]\sum_{i=0}^na_i\sum_{k=0}^n \binom ikx^{k}p^{i-k} \\=&\sum_{i=k}^n\binom ikp^{i-k}a_i \end{align*} \]

对于原式，交换求和顺序：

\[\begin{align*} &\sum_{j=0}^n2^j\sum_{k=0}^j\binom jk(-1)^{j-k}a_k \\=&\sum_{k=0}^na_k\sum_{j=k}^n\binom jk(-1)^{j-k}2^j \end{align*} \]

我们发现，这正好和上面的形式相符。如果取

\[F(x)=\sum_{i=0}^n2^jx^j=\frac{(2x)^{n+1}-1}{2x-1} \]

那么上式的后半部分就是 \(F(x-1)\) 里 \(x^k\) 项的系数。

\[F(x-1)=2^{n+1}\frac{(x-1)^{n+1}-1}{2x-3} \]

二项式定理得到分子的各项系数，除以一次式可以直接递推处理。

于是我们 \(\mathcal O(n)\) 地解决了这个问题。

---

但是上面那个做法太吃注意力了！有没有更简单的方法？

考虑 \(S(n,m)\) 关于 \(n\) 的生成函数 \(\frac 1{m!}(\exp x-1)^m\)

原式相当于

\[\begin{align*} &\sum_{i=0}^n\sum_{j=0}^n 2^j[x^i/i!](\exp x-1)^j \\=&\sum_{i=0}^n\sum_{j=0}^n [x^i/i!](2\exp x-2)^j \\=&\sum_{i=0}^n2^{n+1}[x^i/i!]\frac{(\exp x-1)^{n+1}-1}{2\exp x-3} \end{align*} \]

可以发现得到的结果和上面的方法本质上是一样的。

---

code

#include <algorithm>
#include <iostream>
#include <bitset>

const int MOD = 998244353;
auto gmod = [](auto&& x) { return x - (x >= MOD) * MOD; };
struct mi {
	int w; mi(){}
	mi(int _w): w(_w) {}
	inline int load() { return w; }
	inline bool empty() { return !w; }
};
inline mi operator + (const mi& x, const mi& y) { return gmod(x.w + y.w); }
inline mi operator - (const mi& x, const mi& y) { return gmod(x.w + MOD - y.w); }
inline mi operator - (const mi& x) { return x.w ? MOD - x.w : 0; }
inline mi operator * (const mi& x, const mi& y) { return 1ull * x.w * y.w %MOD; }
inline void operator += (mi& x, const mi& y) { x.w = gmod(x.w + y.w); }
inline void operator -= (mi& x, const mi& y) { x.w = gmod(x.w + MOD - y.w); }
inline void operator *= (mi& x, const mi& y) { x.w = 1ull * x.w * y.w %MOD; }
inline mi fpow(mi x, int k) { mi r = 1; for(; k; k >>= 1, x *= x) if(k & 1) r *= x; return r; }
inline std::istream& operator>> (std::istream& i, mi& m) { i >> m.w; return i; }

namespace # {

const int N = 1e5 + 7;

mi frac[N], ifac[N], inv[N];
inline void init(int n) {
	frac[0] = 1;
	for(int i = 1; i <= n; ++i)
		frac[i] = frac[i-1] * i;
	ifac[n] = fpow(frac[n], MOD-2);
	for(int i = n; i >= 1; --i)
		ifac[i-1] = ifac[i] * i;
	inv[0] = 1;
	for(int i = 1; i <= n; ++i)
		inv[i] = ifac[i] * frac[i-1];
}
inline mi C(int n, int m) { return frac[n] * ifac[m] * ifac[n-m]; }

std::bitset<N> ip;
std::basic_string<int> pr;
mi _pow[N];
void sieve(int n, int power) {
	for(int i = 2; i <= n; ++i) {
		if(!ip[i]) pr += i, _pow[i] = fpow(i, power);
		for(int& p: pr) {
			if(i * p > n) break;
			ip[i * p] = 1, _pow[i * p] = _pow[i] * _pow[p];
			if(i % p == 0) break;
		}
	}
}

int n;
mi tmp[N], poly[N];

inline void main() {
	std::cin >> n;
	init(n+1), sieve(n+1, n+1);	
	mi p2 = fpow(2, n+1);
	for(int i = 0; i <= n + 1; ++i) {
		tmp[i] = p2 * C(n + 1, i);
		if((n + 1 - i) & 1) tmp[i] = -tmp[i];
	}
	tmp[0] -= 1;
	mi inv2 = fpow(2, MOD-2);
	auto calc = [](int k) {
		return k == 0 ? 1 : k == 1 ? n + 1 : (_pow[k] - 1) * inv[k-1];
	};
	mi ans = 0;
	for(int i = n + 1; i >= 1; --i) {
		auto coef = tmp[i] * inv2;
		ans += coef * calc(i-1);
		tmp[i-1] += coef * 3;
	}
	std::cout << ans.load() << "\n";
}

};

int main() {
	#::main();
}