https://gmoj.net/senior/#main/show/6676
考虑用第二类斯特林数把\(i^k\)搞掉：
\(Ans=2^n*f(0)（这一项提出来好，后面不写了）+\sum_{k=1}^m f(k)*\sum_{i=0}^n\binom{n}{i}*i^k\)
\(=\sum_{k=1}^m f(k)*\sum_{i=0}^n\binom{n}{i}*\sum_{j=1}^m \binom{i}{j}*s[k][j]*j!\)
\(=\sum_{k=1}^m f(k)*\sum_{j=1}^m s[k][j]*j!*\sum_{i=0}^n\binom{n}{i}*\binom{i}{j}\)
\(=\sum_{k=1}^m f(k)*\sum_{j=1}^m s[k][j]*j!*\binom{n}{j}*2^{n-j}\)
获得了\(O(m^2)\)的简单做法。
继续搞：
设\(G[j]=j!*\binom{n}{j}*2^{n-j}\)
\(=\sum_{k=1}^m f(k)*\sum_{j=1}^m s[k][j]*G[j]\)
斯特林数有一个容斥的求法：
\(s[n][m]=\frac{1}{m!} \sum_{i=1}^{m(这里填\infty也行)}~~~~ i^n*\binom{m}{i}*(-1)^{m-i}\)
\(原式=\sum_{j=1}^m G(j)*\frac{1}{j!}*\sum_{k=1}^m \sum_{i=1}^m i^k*\binom{j}{i}*(-1)^{j-i}\)
\(=\sum_{j=1}^m G(j)*\frac{1}{j!} * \sum_{i=1}^m\binom{j}{i}*(-1)^{j-i}*\sum_{k=1}^m f(k)*i^k\)
最后一个循环多点求值，再做一遍NTT即可。
Code：
#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i < _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("\n")
using namespace std;
const int mo = 998244353;
ll ksm(ll x, ll y) {
ll s = 1;
for(; y; y /= 2, x = x * x % mo)
if(y & 1) s = s * x % mo;
return s;
}
#define V vector<ll>
#define si size()
#define re resize
namespace ntt {
const int nm = 1 << 18;
ll w[nm], a[nm], b[nm];
int r[nm];
void build() {
for(int i = 1; i < nm; i *= 2) {
w[i] = 1; ll v = ksm(3, (mo - 1) / 2 / i);
ff(j, 1, i) w[i + j] = w[i + j - 1] * v % mo;
}
}
void dft(ll *a, int n, int f) {
ff(i, 0, n) {
r[i] = r[i / 2] / 2 + (i & 1) * (n / 2);
if(i < r[i]) swap(a[i], a[r[i]]);
} ll v;
for(int i = 1; i < n; i *= 2) for(int j = 0; j < n; j += 2 * i) ff(k, 0, i) {
v = a[i + j + k] * w[i + k], a[i + j + k] = (a[j + k] - v) % mo, a[j + k] = (a[j + k] + v) % mo;
}
if(f == -1) {
reverse(a + 1, a + n);
v = ksm(n, mo - 2);
ff(i, 0, n) a[i] = (a[i] + mo) * v % mo;
}
}
V operator * (V p, V q) {
int n0 = p.si + q.si - 1, n = 1;
while(n < n0) n *= 2;
ff(i, 0, n) a[i] = b[i] = 0;
ff(i, 0, p.si) a[i] = p[i];
ff(i, 0, q.si) b[i] = q[i];
dft(a, n, 1); dft(b, n, 1);
ff(i, 0, n) a[i] = a[i] * b[i] % mo;
dft(a, n, -1);
p.re(n0);
ff(i, 0, n0) p[i] = a[i];
return p;
}
void dft(V &p, int f) {
int n = p.si;
ff(i, 0, n) a[i] = p[i];
dft(a, n, f);
ff(i, 0, n) p[i] = a[i];
}
}
using ntt :: operator *;
using ntt :: dft;
V operator - (V a, V b) {
a.re(max(a.si, b.si));
ff(i, 0, a.si) a[i] = (a[i] - b[i]) % mo;
return a;
}
V qni(V a) {
int n0 = a.si;
V b; b.re(1); b[0] = ksm(a[0], mo - 2);
for(int n = 2; n < n0 * 2; n *= 2) {
V c = a; c.re(n); c.re(2 * n);
b.re(2 * n);
dft(c, 1); dft(b, 1);
ff(i, 0, 2 * n) b[i] = (2 * b[i] - c[i] * b[i] % mo * b[i]) % mo;
dft(b, -1); b.re(n);
}
b.re(n0); return b;
}
V qmo(V a, V b) {
int n = a.si - 1, m = b.si - 1;
if(n < m) return a;
V a0 = a, b0 = b;
reverse(a.begin(), a.end());
reverse(b.begin(), b.end());
b.re(a.si);
V c = a * qni(b);
c.re(n - m + 1);
reverse(c.begin(), c.end());
V d = a0 - b0 * c;
d.re(m);
return d;
}
const int N = 1e5 + 5;
int n, m; ll f[N];
ll fac[N], nf[N], h[N];
void build(int n) {
fac[0] = 1; fo(i, 1, n) fac[i] = fac[i - 1] * i % mo;
nf[n] = ksm(fac[n], mo - 2); fd(i, n, 1) nf[i - 1] = nf[i] * i % mo;
}
void build_h() {
h[0] = 1;
fo(i, 1, m) h[i] = h[i - 1] * (n - i + 1) % mo;
fo(i, 0, m) h[i] = h[i] * nf[i] % mo;
}
ll C(int n, int m) {
if(n < m) return 0;
return fac[n] * nf[n - m] % mo * nf[m] % mo;
}
ll g[N];
V t[N * 4];
#define i0 i + i
#define i1 i + i + 1
void dg(int i, int x, int y) {
if(x == y) {
t[i].re(2);
t[i][0] = -x; t[i][1] = 1;
return;
}
int m = x + y >> 1;
dg(i0, x, m); dg(i1, m + 1, y);
t[i] = t[i0] * t[i1];
}
ll p[N];
void fz(V a, int i, int x, int y) {
a = qmo(a, t[i]);
if(y - x + 1 < 128) {
fo(j, x, y) {
ll s = 0, v = 1;
ff(k, 0, a.si) {
s = (s + v * a[k]) % mo;
v = v * j % mo;
}
p[j] = s;
}
return;
}
if(x == y) {
p[x] = a[0];
return;
}
int m = x + y >> 1;
fz(a, i0, x, m); fz(a, i1, m + 1, y);
}
int main() {
ntt :: build();
freopen("number.in", "r", stdin);
freopen("number.out", "w", stdout);
scanf("%d %d", &n, &m);
fo(i, 0, m) scanf("%lld", &f[i]);
build(m);
build_h();
ll ans = f[0] * ksm(2, n) % mo;
fo(j, 1, m) g[j] = fac[j] * ksm(2, n - j) % mo * h[j] % mo;
V a, b; a.clear(); b.clear();
a.re(m + 1); b.re(m + 1);
dg(1, 1, m);
V c; c.re(m + 1);
fo(i, 1, m) c[i] = f[i];
fz(c, 1, 1, m);
fo(i, 1, m) {
a[i] = nf[i] * p[i] % mo;
}
fo(i, 0, m) b[i] = nf[i] * (i % 2 ? -1 : 1);
a = a * b;
fo(i, 1, m) ans = (ans + a[i] * g[i]) % mo;
ans = (ans % mo + mo) % mo;
pp("%lld\n", ans);
}
JZOJ6676. 【2020.06.01省选模拟】查拉图斯特拉如是说 （第二类斯特林数+多项式多点求值）
原文：https://www.cnblogs.com/coldchair/p/13027568.html