斯特林数

简述

这是一篇咕咕咕了三个月的博客，今天来彻底写一写。两类斯特林数在化简式子和本身的组合意义上都对做题有很大的帮助。

第一类斯特林数

组合意义

用 \(n\) 个不同的元素分成 \(k\) 个非空环排列的方案数就是第一类斯特林数，有递推式：

\[ \begin{bmatrix} n \\ k \end{bmatrix} = \begin{bmatrix} n – 1 \\ k – 1 \end{bmatrix} + (n – 1) \begin{bmatrix} n – 1 \\ k \end{bmatrix} \]

解释：考虑独自成环、也可以考虑进入某个环中的某个元素之前。

一些性质

\[ x^{\underline{n}} = \sum_{i = 0}^n \begin{bmatrix} n \\ i \end{bmatrix} (-1)^{n – i} x^i \]

有了上面这个性质，其实做自然数幂和就不难了。设 \(S_k(n) = \sum_{i = 1}^n i^k\)，则有：

\[ \begin{aligned} x^{\underline{n}} &= \sum_{i = 0}^n \begin{bmatrix} n \\ i \end{bmatrix} (-1)^{n – i} x^i \\ x^n &= x^{\underline{n}} – \sum_{i = 0}^{n – 1} \begin{bmatrix} n \\ i \end{bmatrix} (-1)^{n – i} x^i \end{aligned} \]

做和：

\[ \begin{aligned} x^n &= x^{\underline{n}} – \sum_{i = 0}^{n – 1} \begin{bmatrix} n \ i \end{bmatrix} (-1)^{n – i} x^i \\ S_k(n) &= \sum_{i = 1}^n i^k = \sum_{i = 1}^n i^{\underline{k}} – \sum_{j = 0}^{k – 1} \begin{bmatrix} k \ j \end{bmatrix} (-1)^{k – j} i^j \\ &= \frac {(n + 1)^{\underline{k + 1}}} {k + 1} – \sum_{j = 0}^{k – 1} \begin{bmatrix} k \ j \end{bmatrix} (-1)^{k – j} \sum_{i = 1}^n i^j \end{aligned} \]

这个东西用的不多，求自然数幂和用的更多的是下面的第二类斯特林数。

获得方式

考虑用生成函数进行构造：

\[ \sum_{i = 0}^n \begin{bmatrix} n \\ i \end{bmatrix} x^i = \prod_{i = 0}^{n – 1} (x + i) \]

分治 FFT。

// P5408.cpp
#include <bits/stdc++.h>

using namespace std;

const int MAX_N = 6e5 + 20000, mod = 167772161, g = 3;

int n, rev[MAX_N], fac[MAX_N], fac_inv[MAX_N], f[MAX_N];

int quick_pow(int bas, int tim)
{
    int ret = 1;
    while (tim)
    {
        if (tim & 1)
            ret = 1LL * ret * bas % mod;
        bas = 1LL * bas * bas % mod;
        tim >>= 1;
    }
    return ret;
}

const int Gi = quick_pow(g, mod - 2);

void ntt_initialize(int poly_bit)
{
    for (int i = 0; i < (1 << poly_bit); i++)
        rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) << (poly_bit - 1)));
}

void ntt(int *arr, int dft, int poly_siz)
{
    for (int i = 0; i < poly_siz; i++)
        if (i < rev[i])
            swap(arr[i], arr[rev[i]]);
    for (int step = 1; step < poly_siz; step <<= 1)
    {
        int omega_n = quick_pow(dft == 1 ? g : Gi, (mod - 1) / (step << 1));
        for (int j = 0; j < poly_siz; j += (step << 1))
        {
            int omega_nk = 1;
            for (int k = j; k < j + step; k++, omega_nk = 1LL * omega_n * omega_nk % mod)
            {
                int t = 1LL * arr[k + step] * omega_nk % mod;
                arr[k + step] = (0LL + arr[k] + mod - t) % mod, arr[k] = (1LL * arr[k] + t) % mod;
            }
        }
    }
    if (dft == -1)
    {
        int inv = quick_pow(poly_siz, mod - 2);
        for (int i = 0; i < poly_siz; i++)
            arr[i] = 1LL * arr[i] * inv % mod;
    }
}

void offset(int *src, int *dst, int deg, int c)
{
    static int A[MAX_N], B[MAX_N];
    // dst(x) = src(x + c);
    int poly_bit = 0, poly_siz = 1;
    while (poly_siz <= (deg << 1))
        poly_bit++, poly_siz <<= 1;
    for (int i = 0; i < deg; i++)
        A[deg - i - 1] = 1LL * src[i] * fac[i] % mod;
    int pow_acc = 1;
    for (int i = 0; i < deg; i++, pow_acc = 1LL * pow_acc * c % mod)
        B[i] = 1LL * pow_acc * fac_inv[i] % mod;
    for (int i = deg; i < poly_siz; i++)
        A[i] = B[i] = 0;
    ntt_initialize(poly_bit);
    ntt(A, 1, poly_siz), ntt(B, 1, poly_siz);
    for (int i = 0; i < poly_siz; i++)
        A[i] = 1LL * A[i] * B[i] % mod;
    ntt(A, -1, poly_siz);
    for (int i = 0; i < deg; i++)
        dst[i] = 1LL * A[deg - i - 1] * fac_inv[i] % mod;
}

void solve(int deg, int *src)
{
    static int A[MAX_N], B[MAX_N];
    if (deg == 0)
        return (void)(src[0] = 1);
    int len = (deg >> 1);
    solve(len, src);
    offset(src, A, len + 1, len);
    int poly_bit = 0, poly_siz = 1;
    while (poly_siz <= deg)
        poly_bit++, poly_siz <<= 1;
    for (int i = 0; i <= len; i++)
        B[i] = src[i];
    for (int i = len + 1; i < poly_siz; i++)
        A[i] = B[i] = 0;
    ntt_initialize(poly_bit);
    ntt(A, 1, poly_siz), ntt(B, 1, poly_siz);
    for (int i = 0; i < poly_siz; i++)
        A[i] = 1LL * A[i] * B[i] % mod;
    ntt(A, -1, poly_siz);
    if (deg & 1)
        for (int i = 0; i <= deg; i++)
            src[i] = ((i != 0 ? A[i - 1] : 0) + 1LL * (deg - 1) * A[i] % mod) % mod;
    else
        for (int i = 0; i <= deg; i++)
            src[i] = A[i];
}

int main()
{
    scanf("%d", &n);

    for (int i = fac[0] = 1; i <= (n << 1); i++)
        fac[i] = 1LL * fac[i - 1] * i % mod;
    fac_inv[n << 1] = quick_pow(fac[n << 1], mod - 2);
    for (int i = (n << 1) - 1; i >= 0; i--)
        fac_inv[i] = 1LL * fac_inv[i + 1] * (i + 1) % mod;

    solve(n, f);
    for (int i = 0; i <= n; i++)
        printf("%d ", f[i]);
    putchar('\n');
    return 0;
}

第二类斯特林数

组合意义

把 \(n\) 个元素划分成 \(m\) 个子集的方案数就是第二类斯特林数，记作 \( \begin{Bmatrix} n \\ m \end{Bmatrix} \)。

有递推式：

\[ \begin{Bmatrix} n \\ k \end{Bmatrix} = \begin{Bmatrix} n – 1 \\ k – 1 \end{Bmatrix} + k \begin{Bmatrix} n – 1 \\ k \end{Bmatrix} \]

性质

有一个经典的卷积式：

\[ \begin{aligned} \begin{Bmatrix} n \\ k \end{Bmatrix} &= \frac{1}{k!} \sum_{i = 0}^k (-1)^i {k \choose i} (k – i)^n \\ &= \sum_{i = 0}^{k} \frac{(-1)^i}{i!} \frac{(k – i)^n}{(k – i)!} \end{aligned} \]

知道这个之后做一遍多项式乘法就能拿到第二类斯特林数了。

// P5395.cpp
#include <bits/stdc++.h>

using namespace std;

const int MAX_N = 6e5 + 200, mod = 167772161, G = 3;

int n, poly_bit, poly_siz, rev[MAX_N], A[MAX_N], B[MAX_N], res[MAX_N], fac[MAX_N], fac_inv[MAX_N];

int quick_pow(int bas, int tim)
{
    int ret = 1;
    while (tim)
    {
        if (tim & 1)
            ret = 1LL * ret * bas % mod;
        bas = 1LL * bas * bas % mod;
        tim >>= 1;
    }
    return ret;
}

const int Gi = quick_pow(G, mod - 2);

void ntt_initialize()
{
    for (int i = 0; i < poly_siz; i++)
        rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) << (poly_bit - 1)));
}

void ntt(int *arr, int dft)
{
    for (int i = 0; i < poly_siz; i++)
        if (i < rev[i])
            swap(arr[i], arr[rev[i]]);
    for (int step = 1; step < poly_siz; step <<= 1)
    {
        int omega_n = quick_pow(dft == 1 ? G : Gi, (mod - 1) / (step << 1));
        for (int j = 0; j < poly_siz; j += (step << 1))
        {
            int omega_nk = 1;
            for (int k = j; k < j + step; k++, omega_nk = 1LL * omega_nk * omega_n % mod)
            {
                int t = 1LL * omega_nk * arr[k + step] % mod;
                arr[k + step] = (0LL + arr[k] + mod - t) % mod, arr[k] = (0LL + arr[k] + t) % mod;
            }
        }
    }
    if (dft == -1)
    {
        int inv = quick_pow(poly_siz, mod - 2);
        for (int i = 0; i < poly_siz; i++)
            arr[i] = 1LL * arr[i] * inv % mod;
    }
}

int main()
{
    scanf("%d", &n);
    while ((1 << poly_bit) <= (n << 1))
        poly_bit++;
    poly_siz = (1 << poly_bit);
    ntt_initialize();
    for (int i = fac[0] = 1; i <= n; i++)
        fac[i] = 1LL * fac[i - 1] * i % mod;
    fac_inv[n] = quick_pow(fac[n], mod - 2);
    for (int i = n - 1; i >= 0; i--)
        fac_inv[i] = 1LL * fac_inv[i + 1] * (i + 1) % mod;
    for (int i = 0, opt = 1; i <= n; i++)
        A[i] = 1LL * opt * fac_inv[i] % mod, opt = (mod - opt) % mod;
    for (int i = 0; i <= n; i++)
        B[i] = 1LL * quick_pow(i, n) * fac_inv[i] % mod;
    ntt(A, 1), ntt(B, 1);
    for (int i = 0; i < poly_siz; i++)
        res[i] = 1LL * A[i] * B[i] % mod;
    ntt(res, -1);
    for (int i = 0; i <= n; i++)
        printf("%d ", res[i]);
    putchar('\n');
    return 0;
}

还有一个比较重要的东西，就是用第二类斯特林数求自然数幂和。你用二项式反演一下，就可以拿到：

\[ \begin{aligned} n^k &= \sum_{i = 0}^k \begin{Bmatrix} k \\ i \end{Bmatrix} {n \choose i} i! \\ S_k(n) &= \sum_{i = 0}^n \sum_{j = 0}^k \begin{Bmatrix} k \\ j \end{Bmatrix} {i \choose j} j! \\ &= \sum_{j = 0}^k \begin{Bmatrix} k \\ j \end{Bmatrix} j! \sum_{i = 0}^n {i \choose j} \\ &= \sum_{j = 0}^k \begin{Bmatrix} k \\ j \end{Bmatrix} j! {n + 1 \choose j + 1} \end{aligned} \]

一些习题