贝尔数

What · 是什么

贝尔数 B_n（以 Eric Temple Bell 命名，OEIS A000110）是基数为 n 的集合的划分方法总数：把 n 个两两不同的元素分成任意个（可多个）互不区分的非空子集的方案数。序列：B_0=1,B_1=1,B_2=2,B_3=5,B_4=15,B_5=52,B_6=203,…（如 B_3=5：{{1},{2},{3}},{{1},{2,3}},{{2},{1,3}},{{3},{1,2}},{{1,2,3}}）。关键恒等式：B_{n+1}=Σ_{k=0}^n C(n,k)B_k（考虑元素 n+1 与哪 k 个元素同块）；以及 B_n=Σ_{k=0}^n S(n,k)（第二类斯特林数之和）。指数生成函数 EGF：Σ_n B_n x^n/n! = e^{e^x−1}，据此可用多项式 exp O(n log n) 求前 n 项。贝尔三角形（Bell triangle）给出 O(n²) 递推：a_{n,0}=a_{n−1,n−1}，a_{n,m}=a_{n,m−1}+a_{n−1,m−1}，每行首项即贝尔数。模 p 时贝尔数可由 Touchard 同余 B_{n+p} ≡ B_n+B_{n+1} (mod p) 加速。

Why · 为什么

How · 实现

伪代码

# B(n) 贝尔数

def bell_sum_stirling(n, MOD):
  # 用第二类斯特林数行求和: O(n^2)
  S = [[0]*(n+1) for _ in range(n+1)]; S[0][0]=1
  for i in range(1, n+1):
    for j in range(1, i+1):
      S[i][j] = (S[i-1][j-1] + j*S[i-1][j]) % MOD
  return sum(S[n]) % MOD

def bell_triangle(n, MOD):
  # 贝尔三角形 O(n^2)
  row = [1]                       # 第0行
  for i in range(1, n+1):
    nxt = [row[-1]]               # 首项 = 上行末项
    for m in range(1, i+1):
      nxt.append((nxt[m-1] + row[m-1]) % MOD)
    row = nxt
  return row[0]                   # 每行首项是贝尔数

def bell_recurrence(n, MOD):
  # B_{n+1}=sum C(n,k) B_k: 需要组合数, O(n^2) 或 O(n)

Python 实现

MOD = 998244353


def bell_via_stirling(n, mod=MOD):
    """B(n) = sum_k S(n,k), 第二类斯特林数 O(n^2)."""
    row = [0] * (n + 1)
    row[0] = 1
    for i in range(1, n + 1):
        for j in range(i, 0, -1):
            row[j] = (row[j - 1] + j * row[j]) % mod
        row[0] = 0
    return sum(row) % mod


def bell_triangle(n, mod=MOD):
    """贝尔三角形 O(n^2), 返回 [B_0..B_n]."""
    bell = [1]
    row = [1]
    for i in range(1, n + 1):
        nxt = [row[-1]]
        for m in range(1, i + 1):
            nxt.append((nxt[m - 1] + row[m - 1]) % mod)
        row = nxt
        bell.append(row[0])
    return bell


def bell_recurrence(n, mod=MOD):
    """B_{n+1} = sum_k C(n,k) B_k, O(n^2) 预计算组合数."""
    C = [[0] * (n + 1) for _ in range(n + 1)]
    for i in range(n + 1):
        C[i][0] = C[i][i] = 1
        for j in range(1, i):
            C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % mod
    b = [0] * (n + 1)
    b[0] = 1
    for i in range(n):
        s = 0
        for k in range(i + 1):
            s = (s + C[i][k] * b[k]) % mod
        b[i + 1] = s
    return b[n]


def bell_egf_naive(n, mod=MOD):
    """由 EGF e^{e^x-1} 直接按系数递推(形式幂级数 exp) O(n^2):
    exp(A) 的系数 b_k 满足 k*b_k = sum_{j=0}^{k-1} (j+1)*a_{k-1-j}*b_j,
    其中 A(x)=e^x-1, a_j = 1/j!."""
    fact = [1] * (n + 1)
    for i in range(1, n + 1):
        fact[i] = fact[i - 1] * i % mod
    inv = [1] * (n + 1)
    inv[n] = pow(fact[n], mod - 2, mod)
    for i in range(n, 0, -1):
        inv[i - 1] = inv[i] * i % mod
    a = [inv[i] for i in range(n + 1)]      # a_i = 1/i!  (e^x - 1 的非零项)
    b = [0] * (n + 1)
    b[0] = 1
    for k in range(1, n + 1):
        s = 0
        for j in range(k):
            s = (s + (k - j) * a[k - j] * b[j]) % mod   # (k-j)*a_{k-j}
        b[k] = s * pow(k, mod - 2, mod) % mod
    return [b[i] * fact[i] % mod for i in range(n + 1)]  # 还原成 B_n = n! * b_n


if __name__ == "__main__":
    known = {0: 1, 1: 1, 2: 2, 3: 5, 4: 15, 5: 52, 6: 203}
    for n, v in known.items():
        assert bell_via_stirling(n) % MOD == v, n
        assert bell_triangle(n)[n] % MOD == v
        assert bell_recurrence(n) % MOD == v
    egf_list = bell_egf_naive(6)
    for n, v in known.items():
        assert egf_list[n] % MOD == v, (n, egf_list[n])
    assert bell_triangle(6)[4] == 15
    print("bell-number: all tests passed")

易错点

• EGF 系数是 B_n/n!，还原时忘乘 n!
• B_{n+1}=ΣC(n,k)B_k 中求和到 k=n 而非 k=n−1
• 贝尔三角形行内递推 a_{n,m}=a_{n,m−1}+a_{n−1,m−1} 的索引写错
• 把 B_n 与 p(n)（整数分拆）混淆——贝尔数是集合划分
• 模素数用 Touchard 同余时忽略要求 p 为素数

例题

信源

• OI Wiki✓ https://oi-wiki.org/math/combinatorics/bell/
• oiwiki-raw✓ https://raw.githubusercontent.com/OI-wiki/OI-wiki/master/docs/math/combinatorics/bell.md
• oeis✓ https://oeis.org/A000110
• Wikipedia✓ https://en.wikipedia.org/wiki/Bell_number