1. 程式人生 > 其它 >斯特林數入門

斯特林數入門

第二類斯特林數·行

定義:\(\begin{Bmatrix}i\\j\end{Bmatrix}\) 表示將 \(i\) 個元素放入 \(j\) 個不區分的盒子裡的方案數。

遞推:\(\begin{Bmatrix}i\\j\end{Bmatrix} = \begin{Bmatrix}i-1\\j-1\end{Bmatrix} + j\begin{Bmatrix}i-1\\j\end{Bmatrix}\)

二項式反演推導:

首先有 \(m^n = \sum\limits_{i=0}^m\dbinom{m}{i}\begin{Bmatrix}n\\i\end{Bmatrix}i!\)

\(f(i) = i^n,g(i) = \begin{Bmatrix}n\\i\end{Bmatrix}i!\)

則有: \(f(i) = \sum\limits_{j=0}^i\dbinom{i}{j}g(j)\)

那麼有:

\[\begin{aligned}g(i) &= \sum\limits_{j=0}^i(-1)^{i-j}\dbinom{i}{j}f(j)\\&=i!\sum\limits_{j=0}^i\dfrac{(-1)^{i-j}}{(i-j)!} \dfrac{j^n}{j!}\end{aligned} \]

\(\begin{Bmatrix}n\\i\end{Bmatrix} = \sum\limits_{j=0}^i\dfrac{(-1)^{i-j}}{(i-j)!} \dfrac{j^n}{j!}\)

直接 \(ntt\) 捲起來可以做到 \(\mathrm{O(n \log n)}\)

\(\texttt{Code:}\)

//O(n^2)
#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define mp make_pair
#define fi first
#define se second

using namespace std;

inline void read(int &x) {
    x = 0; int f = 0; char c = getchar();
    while (!isdigit(c)) f |= c == '-', c = getchar();
    while (isdigit(c)) x = x * 10 + (c ^ 48), c = getchar();
    x = f ? -x : x;
}

const int cmd = 167772161;
const int N = (1 << 19) + 5;

int fpow(int a, int b) {
    int res = 1;
    for (; b; b >>= 1, a = 1ll * a * a % cmd)
        if (b & 1) res = 1ll * res * a % cmd;
    return res;
}

int add(int a, int b) {a += b; return a < cmd ? a : a - cmd;}

int sub(int a, int b) {a -= b; return a < 0 ? a + cmd : a;}

int n, f[N], g[N], h[N], fac[N], ifac[N];

int C(int n, int m) {
    return 1ll * fac[n] * ifac[m] % cmd * ifac[n - m] % cmd;
}

namespace Poly {
    const int L = (1 << 19) + 5;

    int A[L], B[L], rev[L], yg = 3, invyg = fpow(3, cmd - 2), lby;

    void ntt(int *f, int tp) {
        for (int i = 1; i < lby; i++)
            if (i < rev[i]) swap(f[i], f[rev[i]]);
        for (int len = 2; len <= lby; len <<= 1) {
            int buf0 = fpow(tp ? yg : invyg, (cmd - 1) / len);
            for (int s = 0; s < lby; s += len) {
                int buf = 1;
                for (int i = 0; i < (len >> 1); i++) {
                    int cur = 1ll * buf * f[s + (len >> 1) + i] % cmd;
                    f[s + (len >> 1) + i] = sub(f[s + i], cur);
                    f[s + i] = add(f[s + i], cur);
                    buf = 1ll * buf * buf0 % cmd;
                }
            }
        }
    }

    void NTT(int *h, int *f, int *g) {
        for (int i = 1; i < lby; i++)
            rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? lby >> 1 : 0);
        memcpy(A, f, sizeof(int)*lby);
        memcpy(B, g, sizeof(int)*lby);
        ntt(A, 1); ntt(B, 1);
        for (int i = 0; i < lby; i++) A[i] = 1ll * A[i] * B[i] % cmd;
        ntt(A, 0); int invl = fpow(lby, cmd - 2);
        for (int i = 0; i < lby; i++) h[i] = 1ll * A[i] * invl % cmd;
    }
}

using namespace Poly;

int main() {
    scanf("%d", &n);
    fac[0] = ifac[0] = 1;
    for (int i = 1; i <= n; i++) fac[i] = 1ll * fac[i - 1] * i % cmd, ifac[i] = fpow(fac[i], cmd - 2);
    for (int i = 0; i <= n; i++) f[i] = 1ll * fpow(i, n) * ifac[i] % cmd;
    for (int i = 0; i <= n; i++) g[i] = i & 1 ? cmd - ifac[i] : ifac[i];
    lby = 1; for (; lby <= n + n; lby <<= 1);
    NTT(h, f, g);
    for (int i = 0; i <= n; i++) printf("%d ", h[i]);
    return 0;
}