Description：

$R$ 是圆排列数， $D$ 是错位排列数， $S$ 是第二类斯特林数，求： $$ \sum_{i=1}^n\sum_{j=1}^iS_{i,j}D_jR_j $$ $1\leqslant n\leqslant 10^5$

Solution：

交换一下求和号然后把斯特林数展开： $$ \begin{align} &\sum_{j=1}^nD_jR_j\sum_{i=j}^n\frac1{j!}\sum_{k=0}^j(-1)^{j-k}\binom jki^k\\ =&\sum_{j=1}^nD_jR_j\sum_{k=0}^j\frac{(-1)^{j-k}}{(j-k)!}\frac{\begin{align}\sum_{i=j}^ni^k\end{align}}{k!} \end{align} $$ 看上去非常像一个卷积，但是最后面那个等比数列下标从 $j$ 开始，不能直接卷积。

退回到最初的状态，我们发现原式实际上等于这个： $$ \sum_{i=1}^n\sum_{j=1}^nS_{i,j}D_jR_j $$ 然后再推式子就变成了： $$ \sum_{j=1}^nD_jR_j\sum_{k=0}^j\frac{(-1)^{j-k}}{(j-k)!}\frac{\begin{align}\sum_{i=1}^ni^k\end{align}}{k!} $$ 于是 $NTT$ 就行了。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n;

#define MAXN 400010

int D[MAXN],R[MAXN];

int fac[MAXN],inv[MAXN];

#define P 998244353

int power(int a,int b)

{

	int res = 1;

	for(;b > 0;a = 1ll * a * a % P,b = b >> 1)if(b & 1)res = 1ll * res * a % P;

	return res;

}

int calc(int a,int q,int n)

{

	if(q == 0)return a;if(q == 1)return 1ll * a * n % P;

	return 1ll * a * (power(q,n) + P - 1) % P * power(q + P - 1,P - 2) % P;

}

int ww[MAXN << 1],*g = ww + MAXN;

int rev[MAXN];

void NTT(int f[],int l,int type)

{

	for(int i = 0;i < l;++i)

	{

		if(i < rev[i])swap(f[i],f[rev[i]]);

	}

	for(int i = 2;i <= l;i = i << 1)

	{

		int wn = g[type * l / i];

		for(int j = 0;j < l;j += i)

		{

			int w = 1;

			for(int k = j;k < j + i / 2;++k)

			{

				int u = f[k],t = 1ll * w * f[k + i / 2] % P;

				f[k] = (u + t) % P;

				f[k + i / 2] = (u - t + P) % P;

				w = 1ll * w * wn % P;

			}

		}

	}

	if(type == -1)

	{

		int ni = power(l,P - 2);

		for(int i = 0;i < l;++i)f[i] = 1ll * f[i] * ni % P;

	}

	return;

}

int a[MAXN],b[MAXN],res[MAXN];

int main()

{

	scanf("%d",&n);

	fac[0] = 1;for(int i = 1;i <= n;++i)fac[i] = 1ll * fac[i - 1] * i % P;

	inv[n] = power(fac[n],P - 2);for(int i = n - 1;i >= 0;--i)inv[i] = 1ll * inv[i + 1] * (i + 1) % P;

	R[1] = 1;for(int i = 2;i <= n;++i)R[i] = 1ll * R[i - 1] * (i - 1) % P;

	D[1] = 0;D[0] = 1;for(int i = 2;i <= n;++i)D[i] = 1ll * (i - 1) * (D[i - 1] + D[i - 2]) % P;

	int ans = 0;

	for(int i = 0;i <= n;++i)a[i] = 1ll * power(P - 1,i) % P * inv[i] % P;

	for(int i = 0;i <= n;++i)b[i] = 1ll * calc(i,i,n) * inv[i] % P;

	int l = 0,len = 1;

	for(;len <= n * 2;len = len << 1)++l;

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i < len;++i)g[i] = g[i - len] = 1ll * g[i - 1] * g[1] % P;

	NTT(a,len,1);NTT(b,len,1);

	for(int i = 0;i < len;++i)res[i] = 1ll * a[i] * b[i] % P;

	NTT(res,len,-1);

	for(int j = 1;j <= n;++j)

	{

		ans = (ans + 1ll * res[j] * D[j] % P * R[j] % P) % P;

	}

	cout << ans << endl;

	return 0;

}