求和（第二斯特林数+ntt）

题目：

https://www.luogu.com.cn/problem/P4091

思路：

代码：

#define _CRT_SECURE_NO_WARNINGS
#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<utility>
#include<stack>
#include<queue>
#include<vector>
#include<set>
#include<map>
#include<unordered_map>
using namespace std;
#define  LL long long 
const int N =4e5+100;
const double PI = acos(-1);
LL mod = 998244353,gi=3,gg,g[N],f[N],h[N],rel[N];
int n, tot, bit;
LL quick(LL a, LL b, LL mod)
{
    LL ans = 1;
    while (b)
    {  
        if (b & 1) ans = ans * a % mod;
        b >>= 1;
        a = a * a % mod;
    }
    return ans;
}
void init()
{  
    h[0] = 1;
    for (int i = 1; i <= n; i++)
        h[i] = i * h[i - 1] % mod;
    g[0] = f[0] = 1;
    g[1] = -1;
    f[1] = (n + 1)%mod;
    for (int i = 2; i <= n; i++)
    {
        LL x = quick(h[i], mod - 2, mod), y = quick(i, n + 1, mod)-1;
        g[i] = x;
        if (i % 2) g[i] *= -1;
        x =quick((h[i]*(i - 1))%mod, mod - 2, mod);
        f[i] = x * y % mod;
    }
}
void ntt(LL a[], int op)
{
    for (int i = 0; i < tot; i++)
        if (i < rel[i]) swap(a[i], a[rel[i]]);
    for (int m = 2; m <= tot; m <<= 1)
    {
        LL gk = quick(op == 1 ? gi : gg, (mod - 1) / m, mod);
        for (int i = 0; i < tot; i += m)
        {
            LL  g1 = 1;
            for (int j = 0; j < m / 2; j++)
            {
                LL x = a[i + j], y = a[i + j + m / 2] * g1%mod;
                a[i + j] = (x + y) % mod;
                a[i + j + m / 2] = ((x - y) % mod + mod) % mod;
                g1 = g1 * gk % mod;
            }
      }
    }
    if (op == -1)
    {
        LL gk = quick(tot, mod - 2, mod);
        for (int i = 0; i < tot; i++)
            a[i] = a[i] * gk % mod;
    }
}
int main() {
    cin >> n;
    init();
    while ((1 << bit) < 2 * n + 1) bit++;
    tot = 1 << bit;
    for (int i = 0; i < tot; i++) rel[i] = rel[i / 2] / 2 + ((i & 1) ? tot / 2 : 0);
    gg = quick(gi, mod - 2, mod);
    ntt(f, 1);
    ntt(g, 1);
    for (int i = 0; i < tot; i++)
        f[i] = f[i] * g[i] % mod;
    ntt(f, -1);
    LL ans = 0;
   
    for (int i = 0; i <= n; i++)
    {
        ans = (ans + h[i] * quick(2,i,mod) % mod * f[i] % mod) % mod;
    }
    cout << ans << endl;
    return 0;
}