[洛谷P5408]第一类斯特林数·行
给定 \(n(1\leq n\le 262144)\) ，对于所有的整数 \(i\in [0,n]\)，求出 \({n \brack i} \pmod{167772161}\)
题解
考虑第一类斯特林数的生成函数
\[\sum_{k=0}^{n}{n \brack k}x^k=x^{\overline{n}}
\]
其中生成函数的 \(k\) 次项的系数就是我们要求的第一类斯特林数。
考虑倍增。
\[x^{\overline{2n}}=x^{\overline{n}}(x+n)^{\overline{n}}\\
\]
记\(f(x)=x^{\overline{n}}\)。假设我们已经知道了 \(f(x)\)，如何去求出 \(f(x+n)=(x+n)^{\overline{n}}\)。
令 \(a_i=[x^i]f(x)\)，即 \(a_i\) 是多项式 \(f(x)\) 的第 \(i\) 项前的系数。那么
\[f(x)=\sum_{i=0}^{n}a_ix^i\\
f(x+n)=\sum_{i=0}^{n}a_i(x+n)^i=\sum_{i=0}^{n}a_i\sum_{j=0}^i\binom{i}{j}x^jn^{i-j}\\
=\sum_{j=0}^{n}x^j\sum_{i=j}^{n}\binom{i}{j}a_in^{i-j}=\sum_{j=0}^{n}x^j\sum_{i=j}^{n}a_in^{i-j}\frac{i!}{j!(i-j)!}\\
=\sum_{j=0}^{n}\frac{x^j}{j!}\sum_{i=j}^{n}a_ii!\frac{n^{i-j}}{(i-j)!}\\
\]
记 \(A(i)=a_ii!,B(i)=\frac{n^i}{i!}\)，则
\[f(x+n)=\sum_{j=0}^{n}\frac{x^j}{j!}\sum_{i=j}^{n}A(i)B(i-j)=\sum_{j=0}^{n}\frac{x^j}{j!}\sum_{i=0}^{n-j}A(i+j)B(i)\\
\]
记 \(A'(x)=A(n-x)\) 为 \(A(x)\) 的翻转多项式，则
\[f(x+n)=\sum_{j=0}^{n}\frac{x^j}{j!}\sum_{i=0}^{n-j}A'(n-j-i)B(i)
\]
显然这是一个卷积形式。而模数 \(167772161\) 是一个 NTT 模数，所以我们可以使用 NTT 以 \(O(n\log n)\) 的时间复杂度求出 \(f(x+n)\)，然后再和 \(f(x)\) 进行一次多项式乘法即可求出 \(x^{\overline{2n}}\)。
根据主定理，有
\[T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)
\]
因此时间复杂度为 \(O(n\log n)\)。
Code
#include <bits/stdc++.h>
using namespace std;
#define RG register int
#define LL long long
const LL MOD = 167772161LL;
const int maxn = 262144;
LL inv[maxn + 5], fact[maxn + 5], finv[maxn + 5];
void init() {
inv[1] = fact[0] = fact[1] = finv[0] = finv[1] = 1;
for (int i = 2;i <= maxn;++i) {
inv[i] = ((-(MOD / i) * inv[MOD % i]) % MOD + MOD) % MOD;
fact[i] = fact[i - 1] * i % MOD;
finv[i] = finv[i - 1] * inv[i] % MOD;
}
return;
}
LL ExPow(LL b, LL n, LL MOD) {
LL x = 1;
LL Power = b % MOD;
while (n) {
if (n & 1) x = x * Power % MOD;
Power = Power * Power % MOD;
n >>= 1;
}
return x;
}
int r[maxn + 5];
int L, limit;
const LL P = 167772161, G = 3, Gi = 55924054;
void NTT(LL* A, int type) {
for (int i = 0; i < limit; i++)
if (i < r[i]) swap(A[i], A[r[i]]);
for (int mid = 1; mid < limit; mid <<= 1) {
LL Wn = ExPow(type == 1 ? G : Gi, (P - 1) / (mid << 1), P);
for (int j = 0; j < limit; j += (mid << 1)) {
LL w = 1;
for (int k = 0; k < mid; k++, w = (w * Wn) % P) {
LL x = A[j + k], y = w * A[j + k + mid] % P;
A[j + k] = (x + y) % P;
A[j + k + mid] = (x - y + P) % P;
}
}
}
}
void Convolution(LL* a, int N, LL* b, LL M, LL* c) {
L = 0; limit = 1;
while (limit <= N + M) limit <<= 1, L++;
for (int i = 0; i < limit; i++) {
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
if (i > N) a[i] = 0;
if (i > M) b[i] = 0;
}
NTT(a, 1); NTT(b, 1);
for (int i = 0; i < limit; i++) a[i] = (a[i] * b[i]) % P;
NTT(a, -1);
LL inv = ExPow(limit, P - 2, P);
for (int i = 0; i <= N + M; ++i)
c[i] = a[i] * inv % P;
}
LL a[maxn + 5], b[maxn + 5], c[maxn + 5];
void StirlingI(int n) {
if (n == 1) { c[0] = 0;c[1] = 1;return; }
int m = n >> 1;
StirlingI(m);
LL mi = 1;
for (int i = 0;i <= m;++i) {
a[i] = fact[m - i] * c[m - i] % MOD;
b[i] = mi * finv[i] % MOD;
mi = (mi * m) % MOD;
}
Convolution(a, m, b, m, a);
for (int i = 0;i <= m;++i)
b[i] = finv[i] * a[m - i] % MOD;
Convolution(b, m, c, m, c);
if (n % 2 == 0) return;
c[n] = 0;
for (int i = n;i >= 1;--i)
c[i] = (c[i - 1] + (LL)(n - 1) * c[i] % MOD) % MOD;
c[0] = 0;
return;
}
int main() {
init();
int n;
scanf("%d", &n);
StirlingI(n);
for (int i = 0;i <= n;++i)
printf("%lld ", c[i]);
printf("\n");
return 0;
}
