[洛谷P5408]第一类斯特林数·行

给定 (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;
}