![P4091 [HEOI2016/TJOI2016] 求和](http://pic.xiahunao.cn/yaotu/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)\) 地解决了这个问题。
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();
}