给定 (n(1leq nle 262144)) ，对于所有的整数 (iin [0,n])，求出 ({n rack i} pmod{167772161})
题解
考虑第一类斯特林数的生成函数
[sum_{k=0}^{n}{n rack 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_isum_{j=0}^iinom{i}{j}x^jn^{i-j}\
=sum_{j=0}^{n}x^jsum_{i=j}^{n}inom{i}{j}a_in^{i-j}=sum_{j=0}^{n}x^jsum_{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(nlog n)) 的时间复杂度求出 (f(x+n))，然后再和 (f(x)) 进行一次多项式乘法即可求出 (x^{overline{2n}})。
根据主定理，有
[T(n)=T(frac{n}{2})+O(nlog n)=O(nlog n)
]
因此时间复杂度为 (O(nlog 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("
");
return 0;
}