# 第一类斯特林数的一种 $O(n\log n)$倍增求法

### 快速求出 $S_n^{0\cdots n}$

$$\prod_{i=0}^{n-1}(x+i)=\sum_{i=0}^{n}f_ix^i$$

\begin{aligned} & \prod_{i=n}^{2n-1}(x+i)\\ =& \prod_{i=0}^{n-1}(x+i+n)\\ =& \sum_{i=0}^{n}f_i(x+n)^i \end{aligned}

\begin{aligned} =& \sum_{i=0}^{n}f_i\sum_{j=0}^{i}C_i^jx^jn^{i-j}\\ =& \sum_{j=0}^{n}\sum_{i=j}^{n}f_i\frac{i!}{j!(i-j)!}x^jn^in^{-j}\\ =& \sum_{j=0}^{n}\frac{1}{j!n^j}(\sum_{i=j}^{n}(f_ii!n^i)\frac{1}{(i-j)!})x^j \end{aligned}

• 卷积的下标范围为 $0$到 $n$的闭区间，因此 NTT 时需要获得至少 $2n+1$个点值。
• 求 $\frac{S_{2n}}{S_n}$时，设 $F(x)=\sum_{i=0}f_ii!n^ix^i$，$G(x)=\sum_{i=0}\frac{1}{i!}x^i$。由于公式中卷积的形式与我们平常见到的卷积方向相反，所以我们需要将 $F$翻转，与 $G$相乘，再将结果翻转，才能求出正确的结果。
• 需要保证每次傅里叶变换前数组的高位都是 $0$。

### 代码 (CF960G)

#include <bits/stdc++.h>
#define MOD 998244353
#define G 3
#define GI 332748118LL
#define MOD2 996491788296388609LL
#define MX 262144

using namespace std;

typedef long long ll;

ll qm(const ll& x) {return (x>=MOD) ? (x-MOD) : (x);}

ll qpow(ll x, ll t)
{
ll ans = 1;
while(t)
{
if(t & 1) ans = ans * x % MOD;
x = x * x % MOD;
t >>= 1;
}
return ans;
}

ll inv(ll x)
{
return qpow(x, MOD-2);
}

struct Fourier
{
int n, bit, rev[MX];

void init(int x)
{
n = 1, bit = 0;
while(n < x) n <<=1, ++bit;
for(int i=1; i<n; i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(bit-1));
}

void dft(ll* x, int f)
{
for(int i=0; i<n; i++) if(i < rev[i]) swap(x[i], x[rev[i]]);
for(int w=1,b=1; w<n; w<<=1,b++)
{
ll wx = (f==1 ? qpow(G, (MOD-1)>>b) : qpow(GI, (MOD-1)>>b));
for(int j=0; j<n; j+=(w<<1))
{
ll wn = 1;
for(int i=j; i<j+w; i++)
{
ll a = x[i], b = wn*x[i+w];
x[i] = (a+b) % MOD;
x[i+w] = (a-b+MOD2) % MOD;
wn = wn * wx % MOD;
}
}
}
if(f == -1)
{
ll mul = inv(n);
for(int i=0; i<n; i++) x[i] = x[i] * mul % MOD;
}
}
} NTT;

ll fac[MX], faci[MX];
ll t1[MX], t2[MX], t3[MX], t4[MX];

void get_s(ll* x, int n)
{
if(n == 1)
{
x[0] = 0;
x[1] = 1;
}
else if(n & 1)
{
get_s(x, n-1);
x[n] = 0;
for(int i=n-1; i>=0; i--) x[i+1] = ((n-1) * x[i+1] + x[i]) % MOD;
}
else
{
get_s(x, n/2);
NTT.init(n+1);
for(int i=n/2+1; i<NTT.n; i++) t1[i] = t2[i] = t3[i] = t4[i] = 0;
for(int i=0; i<=n/2; i++)
{
t1[i] = x[i] * qpow(n/2, i) % MOD * fac[i] % MOD;
t2[i] = faci[i];
}
reverse(t1, t1+n/2+1);
NTT.dft(t1, +1);
NTT.dft(t2, +1);
for(int i=0; i<NTT.n; i++) t3[i] = t1[i] * t2[i] % MOD;
NTT.dft(t3, -1);
reverse(t3, t3+n/2+1);
for(int i=0; i<=n/2; i++) t4[i] = t3[i] * inv(qpow(n/2, i)) % MOD * faci[i] % MOD;
NTT.dft(t4, +1);
NTT.dft(x, +1);
for(int i=0; i<NTT.n; i++) x[i] = x[i] * t4[i] % MOD;
NTT.dft(x, -1);
}
}

void init()
{
fac[0] = faci[0] = 1;
for(int i=1; i<MX; i++) fac[i] = fac[i-1] * i % MOD;
faci[MX-1] = inv(fac[MX-1]);
for(int i=MX-1; i>=1; i--) faci[i-1] = faci[i] * i % MOD;
}

ll stiring[MX];
int n, a, b;

int main()
{
scanf("%d%d%d", &n, &a, &b);
if(n==1 && a==1 && b==1) printf("%d\n", 1);
else if(a+b-1>n || a<1 || b<1) printf("%d\n", 0);
else
{
init();
get_s(stiring, n-1);
ll way = stiring[a+b-2];
printf("%lld\n", way * fac[a+b-2] % MOD * faci[a-1] % MOD * faci[b-1] % MOD);
}
return 0;
}


### 1 条评论

#### AABB · 2020年1月2日 4:15 下午

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%