从球盒问题到第二类斯特林数

斯特林数是组合数学中一类特殊的数，有着广泛的应用，本文主要讨论第二类斯特林数的推导，性质与应用。

从球盒问题说起

组合问题最基础的模型就是球盒问题了。球盒问题即为 $n$ 个球放入 $m$ 个盒子的方案数。众所周知，球盒问题有 $2 \times 2 \times 2 = 8$ 种类型，分别为：

1. 球相同，盒子相同，不可以有空盒子。
2. 球相同，盒子相同，可以有空盒子。
3. 球相同，盒子不同，不可以有空盒子。
4. 球相同，盒子不同，可以有空盒子。
5. 球不同，盒子相同，不可以有空盒子。
6. 球不同，盒子相同，可以有空盒子。
7. 球不同，盒子不同，不可以有空盒子。
8. 球不同，盒子不同，可以有空盒子。

大量组合数学问题都源于这 $8$ 种基础的模型，让我们一种一种讨论。

1. 球相同，盒子相同，不可以有空盒子。

这个问题可以转化为整数的 $m$ 拆分数，即把 $n$ 分解成 $m$ 个数相加的方案数，记为 $p_k(n)$ 。 有如下递推式成立：

$$p_k(n) = p_{k - 1}(n - 1) + p_k(n - k)$$
其意义为：分别考虑含有 $1$ 的拆分和不含有 $1$ 的拆分。
边界情况： $p_n(n) = p_1(n) = 1$。

2. 球相同，盒子相同，可以有空盒子。

类似上一种情况，我们也可以用整数拆分数来表示。
我们枚举拆分出数的个数，可以得到下式：

$$p(n, m) = \sum_{i = 1}^{\min(n, m)}p_i(n)$$
特别的，对于 $m \geq n$ 的情况，我们可以用整数拆分数 $p(n)$ 来表示答案，即把 $n$ 分解成若干个数相加的方案数。
$$p(n) = \sum_{i = 1}^{n}p_i(n)$$
对于整数拆分数，我们还可以用生成函数的方法来计算。我们枚举每一个数出现的次数，即可得到整数拆分数的生成函数：
$$\sum_{n = 0}^\infty p(n)x^n = \prod_{i = 1}^\infty(1 + x^i + x^{2i} + x^{3i} + \dots)$$
$$\sum_{n = 0}^\infty p(n)x^n = \prod_{i = 1}^\infty \frac{1}{1 - x^i}$$

3. 球相同，盒子不同，不可以有空盒子。

这是一种比较简单的情况。我们可以用隔板法解决，转化为在 $n$ 个球中间 $n - 1$ 个空隙插入 $m - 1$ 个隔板。方案数为：

$$\binom{n - 1}{m - 1}$$

4. 球相同，盒子不同，可以有空盒子。

类似上一种情况，我们也可以用隔板法，不同的是，这种情况隔板可以相邻。我们有两种思路，一种是预先在每个盒子种放入一个球（即球数增加 $m$ 个）转化为不可以有空盒子的情况，另一种方案是考虑一个 $01$ 序列，$0$ 代表球，$1$ 代表隔板，隔板共 $m - 1$ 个。即在长 $n + m - 1$ 的 $01$ 序列选出 $m - 1$ 个位置放 $1$ ， 其余放 $0$ 的方案数。两种思路都可以得到下式：

$$\binom{n + m - 1}{m - 1}$$

5. 球不同，盒子相同，不可以有空盒子。

这个问题即为 $n$ 个元素划分为 $k$ 个等价集合的方案数，这就是我们今天的重点——第二类斯特林数的定义。
第二类斯特林数通常记作 $S(n, m)$ 或 $\begin{Bmatrix}n \\ m \end{Bmatrix}$ 。
其递推式为：

$$\begin{Bmatrix}n \\ m \end{Bmatrix} = \begin{Bmatrix}n - 1 \\ m - 1 \end{Bmatrix} + m\begin{Bmatrix}n - 1 \\ m \end{Bmatrix}$$
其意义为：考虑第 $n$ 个物品，$n$可以单独构成一个非空集合，此时前 $n-1$ 个物品构成 $m-1$ 个非空的不可辨别的集合，有 $\begin{Bmatrix}n - 1 \\ m - 1 \end{Bmatrix}$ 种方法；也可以前 $n-1$ 种物品构成 $m$ 个非空的不可辨别的集合，第 $n$ 个物品放入任意一个中，这样有 $m\begin{Bmatrix}n - 1 \\ m \end{Bmatrix}$ 种方法。

6. 球不同，盒子相同，可以有空盒子。

对于这一种情况，我们可以枚举几个盒子非空，即为：

$$\sum_{i = 1}^{\min(n, m)} \begin{Bmatrix}n \\ i \end{Bmatrix}$$

7. 球不同，盒子不同，不可以有空盒子。

对于这一种情况，我们发现，在第 $5$ 种情况里的每一个方案，我们把每一个盒子标号，就对应着 $m!$ 种方案。所以这种情况的答案为：

$$\begin{Bmatrix}n \\ m \end{Bmatrix} m!$$

8. 球不同，盒子不同，可以有空盒子。

终于到了最后一种情况，也是最简单的一种。
对于这种情况，我们可以枚举每一个球放在哪个盒子里，就得到了答案：

$$m^n$$
当然我们也可以枚举哪一些盒子非空，就得到了一个重要的等式：
$$m^n = \sum_{i = 0}^m \binom{m}{i}\begin{Bmatrix}n \\ i \end{Bmatrix} i!$$
关于这个等式的应用将在之后介绍。

以上我们解决了球盒问题的所有 $8$ 种情况，现在看来，是不是很 $naive$ 呢？

第二类斯特林数

第二类斯特林数是 $n$ 个元素划分为 $k$ 个等价集合的方案数。也就是之前介绍的球盒问题的第 $5$ 种类型。
第二类斯特林数除了解决组合问题，还有一些神奇的性质。
第二类斯特林数可以使用下降幂来表示数的幂，即满足下式：

$$x^n = \sum_{k = 0}^n\begin{Bmatrix}n \\ k \end{Bmatrix} x^{\underline{k}}$$
其中 $x^{\underline{k}} = x(x-1)(x-2)\dots(x-k+1) = \frac{x!}{(x-k)!}$ 。
这个式子可以使用数学归纳法证明，但是我们惊奇的发现，他所表示的意义即为 球不同，盒子不同，可以有空盒子 中我们发现的等式！
把下降幂用组合数替代更直观一些：

$$x^n = \sum_{k = 0}^n\begin{Bmatrix}n \\ k \end{Bmatrix} \frac{x!}{k!(x-k)!} k! = \sum_{k = 0}^n\begin{Bmatrix}n \\ k \end{Bmatrix} \binom{x}{k}k!$$

我们还可以利用这个性质，再次得到第二类斯特林数的递推公式：

转自riteme的博客

$$\begin{aligned} n^{p-1} & = \sum_{k = 0}^{p-1} \begin{Bmatrix}p-1 \\ k \end{Bmatrix}n^{\underline{k}} \\ n^p = n \cdot n^{p-1} & = n\sum_{k = 0}^{p-1} \begin{Bmatrix}p-1 \\ k \end{Bmatrix}n^{\underline{k}} \\ & = \sum_{k = 0}^{p-1} \begin{Bmatrix}p-1 \\ k \end{Bmatrix} \cdot n \cdot n^{\underline{k}} \\ & = \sum_{k = 0}^{p-1} \begin{Bmatrix}p-1 \\ k \end{Bmatrix} \cdot (n - k + k) \cdot n^{\underline{k}} \\ & = \sum_{k = 0}^{p-1} \begin{Bmatrix}p-1 \\ k \end{Bmatrix} n^{\underline{k + 1}} + \sum_{k = 0}^{p-1} k \begin{Bmatrix}p-1 \\ k \end{Bmatrix} n^{\underline{k}} \\ & = \sum_{k = 1}^{p} \begin{Bmatrix}p-1 \\ k - 1\end{Bmatrix} n^{\underline{k}} + \sum_{k = 0}^{p-1} k \begin{Bmatrix}p-1 \\ k \end{Bmatrix} n^{\underline{k}} \end{aligned}$$
我们就再次得到了第二类斯特林数的递推式：
$$\begin{Bmatrix}p \\ k \end{Bmatrix} = k\begin{Bmatrix}p-1 \\ k \end{Bmatrix} + \begin{Bmatrix}p-1 \\ k-1 \end{Bmatrix}$$

现在，你应该对第二类斯特林数有了深入的理解了。

斯特林反演

根据第二类斯特林数的递推式，我们已经可以以 $O(n^2)$ 的复杂度求第二类斯特林数了。那么对于 $n,m$ 比较大的情况，我们应该怎么办呢？这里介绍一个技巧——斯特林反演。

斯特林反演的公式如下：

$$\begin{Bmatrix}n \\ j \end{Bmatrix} = \frac{1}{j!} \sum_{k = 0}^{j} (-1)^{j - k} \binom{j}{k} k^n$$
看起来十分玄学，我们需要想办法证明这个式子。
我们首先对于等号两边乘以阶乘：
$$\begin{Bmatrix}n \\ j \end{Bmatrix} j! = \sum_{k = 0}^{j} (-1)^{j - k} \binom{j}{k} k^n$$
考虑左式的组合意义： $\begin{Bmatrix}n \\ m \end{Bmatrix}m!$ 表示球盒问题中 球不同，盒子不同，不可以有空盒子 的方案数。
我们考虑用另一种思路求解这个组合问题：
考虑容斥原理，没有空盒子的方案数 = 总方案数 - 至少一个空盒子的方案数 + 至少两个个空盒子的方案数 - 至少三个空盒子的方案数……
即
$$\sum_{k = 0}^m (-1)^{m - k} \binom{m}{k} k^n$$
恰好与等号右边相同！
于是我们证明了斯特林反演。
我们把斯特林反演的式子展开：
$$\begin{Bmatrix}n \\ j \end{Bmatrix} =\frac{1}{j!} \sum_{k = 0}^{j} (-1)^{j - k} \frac{j!}{k!(j-k)!} k^n$$
$$\begin{Bmatrix}n \\ j \end{Bmatrix} =\sum_{k = 0}^{j} \frac{(-1)^{j - k}}{(j-k)!} \frac{k^n}{k!}$$
恰好是一个卷积的形式。于是我们可以用 FFT（NTT）以 $O(n \log_2 n)$ 的复杂度算出一行的第二类斯特林数。

例题

BZOJ 2159 Crash 的文明世界

题意：

给定一棵 $n$ 个点的树和一个常数 $k$ , 对于每个 $i$ , 求

$$S(i) = \sum_{j = 1}^n \mathrm{dist}(i, j)^k$$
$n \leq 50000, k \leq 150$ 。

题解：

常见套路：先大力把幂转为下降幂，再转为组合数：

$$\begin{align} S(i) & = \sum_{j = 1}^n \sum_{l = 0}^k \begin{Bmatrix}k \\ l \end{Bmatrix} \mathrm{dist}(i, j)^{\underline{l}}\\ & = \sum_{l = 0}^k \begin{Bmatrix}k \\ l \end{Bmatrix} \sum_{j = 1}^n \mathrm{dist}(i, j)^{\underline{l}}\\ & = \sum_{l = 0}^k \begin{Bmatrix}k \\ l \end{Bmatrix} l! \sum_{j = 1}^n \binom{\mathrm{dist}(i, j)}{l}\\ \end{align}$$
然后根据组合数的递推公式
$$\binom{n}{m} = \binom{n - 1}{m - 1} + \binom{n - 1}{m}$$
就可以 $O(k)$ 从一个点转移到另一个点了。
总时间复杂度 $O(nk)$ 。

My Code：

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <queue>
#include <vector>
#define inf 0x3f3f3f3f

using namespace std;
typedef long long ll;
const ll mod = 10007;

struct edge{
    int to, nxt;
}e[100005];

int h[50005], cnt;

void addedge(int x, int y){
    cnt++; e[cnt].to = y; e[cnt].nxt = h[x]; h[x] = cnt;
    cnt++; e[cnt].to = x; e[cnt].nxt = h[y]; h[y] = cnt;    
}
ll f[50005][151];
int n, k, L;

void dfs(int x, int fa){
    f[x][0] = 1;
    for(int i = h[x]; i; i = e[i].nxt){
        if(e[i].to == fa) continue;
        dfs(e[i].to, x);
        f[x][0] = (f[x][0] + f[e[i].to][0]) % mod;
        for(int j = 1; j <= k; j ++){
            f[x][j] = (f[x][j] + f[e[i].to][j] + f[e[i].to][j - 1]) % mod;
        }
    }
}

void dfs2(int x, int fa){
    for(int i = h[x]; i; i = e[i].nxt){
        if(e[i].to == fa) continue;
        for(int j = k; j >= 2; j --){
            f[e[i].to][j] = (f[x][j] + f[x][j - 1] - 2 * f[e[i].to][j - 1] - f[e[i].to][j - 2] + mod * 3) % mod;
        }
        f[e[i].to][1] = (f[x][1] + f[x][0] - 2 * f[e[i].to][0] + mod + mod) % mod;
        f[e[i].to][0] = f[x][0];
        dfs2(e[i].to, x);
    }
}
ll S[205][205];
ll fac[205];
int main(){
    scanf("%d%d%d", &n, &k, &L);
    int now, a, b, q;
    scanf("%d%d%d%d", &now, &a, &b, &q);
    for(int i = 1; i < n; i ++){
        now = (now * a + b) % q;
        int tmp=((i<L)?i:L);
        int u=i-now%tmp; int v=i+1;
        addedge(u, v);
    }
    dfs(1, 0);
    dfs2(1, 0);
    S[0][0] = 1;
    for(int i = 1; i <= k; i ++){
        S[i][0] = 0;
        for(int j = 1; j <= i; j ++){
            S[i][j] = (S[i - 1][j - 1] + j * S[i - 1][j]) % mod;
        }
    }
    fac[0] =1;
    for(int i = 1; i <= k; i ++){
        fac[i] = fac[i - 1] * i % mod;
    }
    for(int i = 1; i <= n; i ++){
        ll ans = 0;
        for(int j = 0; j <= k; j ++){
            ans = (ans + S[k][j] * fac[j] % mod * f[i][j])% mod;
        }
        printf("%lld\n", ans);
    }
    return 0;
}

HNOI2015 省队集训 原题的价值

题意：

定义一个无向图的权值为所有节点的度数的 $k$ 次方之和. 求所有 $n$个点的简单无向图（无重边无自环）的权值之和, 对 $998244353$ 取模。
$n \leq 10^9,k \leq 2 \times 10^5$ 。

题解：

由于每个点是等价的，我们不妨考虑一个点对答案的贡献，最后乘 $n$ 即为答案。
我们枚举这个点连出去的 $n - 1$ 条边是否存在，其他的部分（一个 $n - 1$ 个点的无向图）与这个点的度数无关，所以再乘上 $n - 1$ 个点的简单无向图的数量。所以答案为：

$$\left ( \sum_{i = 0}^{n - 1}\binom{n - 1}{i}i^k \right ) 2^{\frac{(n - 1)(n - 2)}{2}}n$$
问题转化为求
$$\sum_{i = 0}^n \binom{n}{i}i^k$$
$$\begin{align} \sum_{i = 0}^n \binom{n}{i}i^k &= \sum_{i = 0}^n \binom{n}{i} \sum_{j = 0}^k \begin{Bmatrix}k \\ j \end{Bmatrix} i^{\underline{j}} \\ &= \sum_{i = 0}^n \binom{n}{i} \sum_{j = 0}^k \begin{Bmatrix}k \\ j \end{Bmatrix} \binom{i}{j}j! \\ &= \sum_{j = 0}^k \begin{Bmatrix}k \\ j \end{Bmatrix} j! \sum_{i = 0}^n \binom{n}{i} \binom{i}{j} \\ &= \sum_{j = 0}^k \begin{Bmatrix}k \\ j \end{Bmatrix} j! \sum_{i = 0}^n \frac{n!}{i!(n - i)!} \frac{i!}{j!(i - j)!} \\ &= \sum_{j = 0}^k \begin{Bmatrix}k \\ j \end{Bmatrix} j! \sum_{i = 0}^n \frac{n!}{j!(n - j)!}\frac{(n - j)!}{(n - i)!(i - j)!} \\ &= \sum_{j = 0}^k \begin{Bmatrix}k \\ j \end{Bmatrix} j! \binom{n}{j} \sum_{i =j}^n \binom{n - j}{i - j} \\ &= \sum_{j = 0}^k \begin{Bmatrix}k \\ j \end{Bmatrix} n^{\underline{j}} 2^{n - j} \\ \end{align}$$
这样已经可以 $O(k^2)$ 的递推斯特林数求解了。 弱化版提交地址
观察到我们只需要一行斯特林数，我们可以用斯特林反演 + NTT 在 $O(n \log_2 n)$ 的时间求出答案。

My Code：

#include <bits/stdc++.h>
#define MAXN 524288
using namespace std;
typedef long long ll;
const ll mod = 998244353;

ll qpow(ll a, ll b){
    ll ret = 1;
    for(; b; b >>= 1, a = a * a % mod){
        if(b & 1) ret = ret * a % mod;
    }
    return ret;
}

ll inv3, pw3[MAXN + 1], ipw3[MAXN + 1];
ll fac[MAXN], ifac[MAXN];
int N;
void pre(){
    for(int i = 1; i <= N; i <<= 1) pw3[i] = qpow(3, (mod - 1) / i);
    inv3 = qpow(3, mod - 2);
    for(int i = 1; i <= N; i <<= 1) ipw3[i] = qpow(inv3, (mod - 1) / i);
    fac[0] = 1;
    for(int i = 1; i < N; i ++) fac[i] = fac[i - 1] * i % mod;
    ifac[N - 1] = qpow(fac[N - 1], mod - 2);
    for(int i = N - 1; i >= 1; i --) ifac[i - 1] = ifac[i] * i % mod;
}

void bit_reverse(ll *r, int N){
    for(int i = 0, j = 0; i < N; i ++){
        if(i < j) swap(r[i], r[j]);
        for(int l = N >> 1; (j ^= l) < l; l >>= 1);
    }
}

void NTT(ll *r, int N, int flag){
    bit_reverse(r, N);
    for(int i = 2; i <= N; i <<= 1){
        int m = i >> 1;
        for(int j = 0; j < N; j += i){
            ll w = 1, wn = pw3[i];
            if(flag == -1) wn = ipw3[i];
            for(int k = 0; k < m; k ++){
                ll z = w * r[j + k + m] % mod;
                r[j + k + m] = (r[j + k] - z + mod) % mod;
                r[j + k] = (r[j + k] + z) % mod;
                w = w * wn % mod;
            }
        }
    }
    if(flag == -1){
        ll invn = qpow(N, mod - 2);
        for(int i = 0; i < N; i ++){
            r[i] = r[i] * invn % mod;
        }
    }
}

ll a[MAXN], b[MAXN];
int n, k;
int main(){
    scanf("%d%d", &n, &k);
    N = 1;
    while(N < 2 * k) N <<= 1;
    pre();
    for(int i = 0; i <= k; i ++){
        if(i & 1) a[i] = mod - 1;
        else a[i] = 1;
        a[i] = a[i] * ifac[i] % mod;
        b[i] = qpow(i, k) * ifac[i] % mod;
    }
    NTT(a, N, 1);
    NTT(b, N, 1);
    for(int i = 0; i < N; i ++) a[i] = a[i] * b[i] % mod;
    NTT(a, N, -1);
    n--;
    ll ans = 0, tmp = qpow(2, n) % mod, inv2 = qpow(2, mod - 2);
    for(int i = 0; i <= k; i ++){
        ans = (ans + tmp * a[i]) % mod;
        tmp = tmp * (n - i) % mod * inv2 % mod;
    }
    n++;
    ans = ans * qpow(2, 1ll * (n - 1) * (n - 2) / 2) % mod * n % mod;
    printf("%lld\n", ans);
    return 0;
}


---

参考资料：

斯特灵数 - 维基百科
差分序列与Stirling数 - riteme的博客