Definition

Stirling numbers are fundamental combinatorial objects categorized into two distinct types:

• First kind: Counts the number of ways to partition a set of $n$ elements into $m$ circular permutations.
• Second kind: Counts the number of ways to partition a set of $n$ elemants into $m$ non-empty unordered subsets.

Recurrence Relations

Both types follow recurrence relations derived from combinatorial analysis:

• First kind: $s(n, m) = s(n-1, m-1) + (n-1) \cdot s(n-1, m)$
• Second kind: $S(n, m) = S(n-1, m-1) + m \cdot S(n-1, m)$

Fast Computation via Polynomial Methods

Standard recurrence computation has $O(nm)$ complexity. Using generating functions and polynomial convolution (via Number Theoretic Transofrm), entire rows or columns can be computed in $O(n \log n)$ time.

Second Kind: Row Computation

The identity $S(n, m) = \frac{1}{m!} \sum_{i=0}^{m} (-1)^{m-i} \binom{m}{i} i^n$ trensforms the problem into polynomial convolution. Define sequences $a_i = \frac{i^n}{i!}$ and $b_i = \frac{(-1)^i}{i!}$, where the result is $a * b$.

#include <cstdio>
#include <algorithm>
using namespace std;

typedef long long ll;

const int MAXN = 1000005;
const ll MOD = 167772161;
const int ROOT = 3;

int set_size;
int rev_arr[MAXN];
ll inv_fact[MAXN];
ll polyA[MAXN], polyB[MAXN], root_powers[MAXN];

void compute_bit_reversal(int n) {
    for (int i = 1; i < n; i++) {
        rev_arr[i] = (rev_arr[i >> 1] >> 1) | ((i & 1) ? (n >> 1) : 0);
    }
}

ll modular_pow(ll base, ll exp, ll mod) {
    ll result = 1;
    while (exp) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

void ntt_transform(ll *arr, int n) {
    compute_bit_reversal(n);
    for (int i = 1; i < n; i++) {
        if (rev_arr[i] > i) swap(arr[i], arr[rev_arr[i]]);
    }
    for (int len = 2, half = 1; len <= n; half = len, len <<= 1) {
        ll wn = modular_pow(ROOT, (MOD - 1) / len, MOD);
        root_powers[0] = 1;
        for (int i = 1; i < len; i++) {
            root_powers[i] = root_powers[i - 1] * wn % MOD;
        }
        for (int l = 0, r = len - 1; r < n; l += len, r += len) {
            int w = 0;
            for (int p = l; p < l + half; p++, w++) {
                ll x = arr[p];
                ll y = root_powers[w] * arr[p + half] % MOD;
                arr[p] = (x + y) % MOD;
                arr[p + half] = (x - y + MOD) % MOD;
            }
        }
    }
}

void inverse_ntt(ll *arr, int n) {
    ntt_transform(arr, n);
    reverse(arr + 1, arr + n);
    ll inv_n = modular_pow(n, MOD - 2, MOD);
    for (int i = 0; i < n; i++) {
        arr[i] = arr[i] * inv_n % MOD;
    }
}

int main() {
    scanf("%d", &set_size);
    inv_fact[0] = inv_fact[1] = 1;
    for (int i = 2; i <= set_size; i++) {
        inv_fact[i] = (MOD - MOD / i) * inv_fact[MOD % i] % MOD;
    }
    for (int i = 1; i <= set_size; i++) {
        inv_fact[i] = inv_fact[i] * inv_fact[i - 1] % MOD;
    }
    int size = set_size + 1;
    for (int i = 0; i < size; i++) {
        polyA[i] = modular_pow(i, set_size, MOD) * inv_fact[i] % MOD;
        polyB[i] = (i & 1) ? (MOD - inv_fact[i]) : inv_fact[i];
    }
    int transform_size = 1;
    while (transform_size < (size << 1)) transform_size <<= 1;
    ntt_transform(polyA, transform_size);
    ntt_transform(polyB, transform_size);
    for (int i = 0; i < transform_size; i++) {
        polyA[i] = polyA[i] * polyB[i] % MOD;
    }
    inverse_ntt(polyA, transform_size);
    for (int i = 0; i < size; i++) {
        printf("%lld ", polyA[i]);
    }
    return 0;
}

First Kind: Row Computation

Using the generating function $x^{\overline{n}} = \sum_{i=0}^{n} {n \brack i} x^i$ and identity $x^{\overline{2n}} = x^{\overline{n}} (x + n)^{\overline{n}}$, we apply recursive doubling with polynomial shifting to compute the row efficiently.

#include <cstdio>
#include <algorithm>
using namespace std;

typedef long long ll;

const int MAXN = 1 << 19;
const ll MOD = 167772161;
const int ROOT = 3;

int n;
int rev_arr[MAXN];
ll fact[MAXN], inv_fact[MAXN];
ll polyA[MAXN], polyB[MAXN], temp_poly[MAXN];

void compute_bit_reversal(int n) {
    for (int i = 1; i < n; i++) {
        rev_arr[i] = (rev_arr[i >> 1] >> 1) | ((i & 1) ? (n >> 1) : 0);
    }
}

ll modular_pow(ll base, ll exp, ll mod) {
    ll result = 1;
    while (exp) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

void ntt(ll *arr, int n) {
    compute_bit_reversal(n);
    for (int i = 1; i < n; i++) {
        if (rev_arr[i] > i) swap(arr[i], arr[rev_arr[i]]);
    }
    for (int len = 2, half = 1; len <= n; half = len, len <<= 1) {
        ll wn = modular_pow(ROOT, (MOD - 1) / len, MOD);
        ll w = 1;
        for (int i = 0; i < len; i++) {
            if (i) w = w * wn % MOD;
            temp_poly[i] = w;
        }
        for (int i = 0; i < n; i += len) {
            for (int j = 0; j < half; j++) {
                ll x = arr[i + j];
                ll y = temp_poly[j] * arr[i + j + half] % MOD;
                arr[i + j] = (x + y) % MOD;
                arr[i + j + half] = (x - y + MOD) % MOD;
            }
        }
    }
}

void inverse_ntt(ll *arr, int n) {
    ntt(arr, n);
    reverse(arr + 1, arr + n);
    ll inv_n = modular_pow(n, MOD - 2, MOD);
    for (int i = 0; i < n; i++) {
        arr[i] = arr[i] * inv_n % MOD;
    }
}

void multiply_poly(ll *a, ll *b, int len, int size) {
    int k = 1;
    while (k < (len << 1)) k <<= 1;
    ntt(a, k);
    ntt(b, k);
    for (int i = 0; i < k; i++) {
        a[i] = a[i] * b[i] % MOD;
    }
    inverse_ntt(a, k);
}

void shift_poly(ll *src, ll *dst, int len, int shift) {
    for (int i = 0; i < len; i++) {
        temp_poly[len - i - 1] = src[i] * fact[i] % MOD;
    }
    ll coef = 1;
    for (int i = 0; i < len; i++, coef = coef * shift % MOD) {
        dst[i] = coef * inv_fact[i] % MOD;
    }
    multiply_poly(temp_poly, dst, len, len);
    for (int i = 0; i < len; i++) {
        dst[len - i - 1] = temp_poly[i] * inv_fact[len - i - 1] % MOD;
    }
}

void compute_stirling_row(ll *stirling, ll *temp, int n) {
    if (n == 1) {
        stirling[0] = 0;
        stirling[1] = 1;
        return;
    }
    if (n & 1) {
        compute_stirling_row(stirling, temp, n - 1);
        stirling[n] = 0;
        for (int i = n; i > 0; i--) {
            stirling[i] = (stirling[i - 1] + (n - 1) * stirling[i]) % MOD;
        }
        stirling[0] = stirling[0] * (n - 1) % MOD;
        return;
    }
    compute_stirling_row(stirling, temp, n >> 1);
    shift_poly(stirling, temp, (n >> 1) + 1, n >> 1);
    multiply_poly(stirling, temp, (n >> 1) + 1, n + 1);
}

int main() {
    scanf("%d", &n);
    fact[0] = inv_fact[0] = 1;
    for (int i = 1; i <= n; i++) {
        fact[i] = fact[i - 1] * i % MOD;
        inv_fact[i] = (MOD - MOD / i) * inv_fact[MOD % i] % MOD;
    }
    for (int i = 1; i <= n; i++) {
        inv_fact[i] = inv_fact[i] * inv_fact[i - 1] % MOD;
    }
    compute_stirling_row(polyA, polyB, n);
    for (int i = 0; i <= n; i++) {
        printf("%lld ", polyA[i]);
    }
    return 0;
}