ベル数を求めるアルゴリズム：第2種スターリング数の和を効率的に計算する
ベル数は「写像12相」と呼ばれるものを通して学ぶと、他の数え上げ問題との関わりが分かり全体像がスッキリします。ぜひ確認してみることをおすすめします。
「写像12相」で典型的な数え上げ問題のパターン総整理
ベル数とは
ベル数とは「区別できる \(n\) 個の要素を区別できない \(k\) 個以下のブロックに分割する方法の数」のことです。
これは、第2種スターリング数を用いて表すことが可能です。
第2種スターリング数とは、区別できる \(n\) 個の要素を区別できない \(k\) 個のブロックに分割する方法の数のことで、\(S(n,k)\) などと表すことが多いです。
ベル数は、第2種スターリング数の和を用いて以下のように表せます。
$$B(n,k) =\sum_{i=0}^{k}S(n,i)$$
\(n=k\) の時に限った \(B(n) =\sum_{i=0}^{n}S(n,i)\) をベル数と呼ぶことも多いです。
アルゴリズム
いくつか方針が考えられますが、ここでは2つの方法を紹介します。
- 単純に第2種スターリング数を求めてから和を計算して \(B(n,k) =\sum_{i=0}^{k}S(n,i)\) を求める
- \(B(n,k) = \sum_{i = 0}^{k} \frac{i^n}{i!} \sum_{j = 0}^{k-i} \frac{(-1)^j}{j!}\) となることを利用して計算する
第2種スターリング数の和を求めるアルゴリズム：
- 必要となる第2種スターリング数を動的計画法などで求めておく
- \(B(n,k) =\sum_{i=0}^{k}S(n,i)\) を求める
※ 第2種スターリング数は、有名公式 \(S(n,k) = k S(n-1,k) + S(n-1, k-1)\) を利用して、動的計画法により求めることができます。
\(\sum_{i = 0}^{k} \frac{i^n}{i!} \sum_{j = 0}^{k-i} \frac{(-1)^j}{j!}\) を求めるアルゴリズム：
- \(\sum_{j = 0}^{k-i} \frac{(-1)^j}{j!}\) を前計算しておく
- \(\sum_{i = 0}^{k} \frac{i^n}{i!} \sum_{j = 0}^{k-i} \frac{(-1)^j}{j!}\) を計算する
※ \(B(n,k) =\sum_{j=0}^{k}S(n,j)\) を以下のように式変形できることを利用しています。
\begin{align} B(n,k) &=\sum_{j=0}^{k}S(n,j) \\ &= \sum_{j = 0}^{k} \frac{1}{j!} \sum_{i = 0}^{j} (-1)^{j-i} {}{j}{\rm C}{i} i^{n} \\ &= \sum_{j = 0}^{k} \sum_{i = 0}^{j} \frac{(-1)^{j-i}}{i!(j-i)!} i^{n} \\ &= \sum_{i = 0}^{k} \sum_{j = i}^{k} \frac{(-1)^{j-i}}{i!(j-i)!} i^{n} \\ &= \sum_{i = 0}^{k} \sum_{j = 0}^{k-i} \frac{(-1)^{j}}{i! j!} i^{n}\\ &= \sum_{i = 0}^{k} \frac{i^n}{i!} \sum_{j = 0}^{k-i} \frac{(-1)^j}{j!} \end{align}
計算量
- 第2種スターリング数の和を求めるアルゴリズム：\(O(kn)\)
- \(\sum_{i = 0}^{k} \frac{i^n}{i!} \sum_{j = 0}^{k-i} \frac{(-1)^j}{j!}\) を求めるアルゴリズム：\(O(k \log n)\)
第2種スターリング数の和を求めるアルゴリズムでは、\(S(n,k)\) を求めるため動的計画法により \(n \times k\) のテーブルを更新することになります。1つのマスを計算するのに \(O(1)\) かかるので、全体では \(O(nk)\) だけかかってしまいます。
一度スターリング数を求めてしまえば、後はそれを \(k\) 個足すだけなので、この部分は全体の計算量に影響しません。
\(\sum_{i = 0}^{k} \frac{i^n}{i!} \sum_{j = 0}^{k-i} \frac{(-1)^j}{j!}\) を求めるアルゴリズムでは、\(\sum_{j = 0}^{k-i} \frac{(-1)^j}{j!}\) の部分については前計算を行うことで \(O(1)\) 程度で計算できます。\(i^n\) を求めるのには\(O(\log n)\) だけかかるので、\(k\) 回分計算すると全体では \(O(k \log n)\) になります。
以上から分かりますが、実はベル数を求めるのに必要な計算量は第2種スターリング数を求めるための計算量と大きく変わりません。
C++での実装例
AOJ DPL_5_G Balls and Boxes 7 を解く実装例です。ベル数を求めます。
第2種スターリング数の和を求めるアルゴリズム
#include <bits/stdc++.h>
using namespace std;
const int MOD = 1000000007; // 素数
int main() {
long long N, K;
cin >> N >> K;
vector<vector<long long>> S(N + 1, vector<long long>(K + 1, 0));
S[0][0] = 1;
for (long long n = 1; n < N + 1; n++) {
for (long long k = 1; k < K + 1; k++) {
S[n][k] = (k * S[n - 1][k] + S[n - 1][k - 1]) % MOD;
}
}
long long ans = 0;
for (int k = 0; k <= K; k++) {
ans = (ans + S[N][k]) % MOD;
}
cout << ans << endl;
return 0;
}
\(\sum_{i = 0}^{k} \frac{i^n}{i!} \sum_{j = 0}^{k-i} \frac{(-1)^j}{j!}\) を求めるアルゴリズム
#include <bits/stdc++.h>
using namespace std;
const int MOD = 1000000007; // 素数
vector<long long> pre_sum, fact_inv, inv;
/* init_bell :前処理
計算量:O(n)
*/
void init_bell(int SIZE) {
// 階乗の逆元を計算
fact_inv.resize(SIZE + 5);
inv.resize(SIZE + 5);
fact_inv[0] = fact_inv[1] = 1;
inv[1] = 1;
for (int i = 2; i < SIZE + 5; i++) {
inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD;
fact_inv[i] = fact_inv[i - 1] * inv[i] % MOD;
}
// 目的の値を前計算
pre_sum.resize(SIZE + 5);
pre_sum[0] = 1;
for (int i = 1; i < SIZE + 5; i++) {
if (i % 2 == 1) {
pre_sum[i] = (pre_sum[i - 1] + (MOD - fact_inv[i])) % MOD;
} else {
pre_sum[i] = (pre_sum[i - 1] + fact_inv[i]) % MOD;
}
}
}
// べき乗
long long pow(long long x, long long n) {
long long ret = 1;
while (n > 0) {
if (n & 1) ret = ret * x % MOD;
x = x * x % MOD;
n >>= 1;
}
return ret;
}
/* ベル数を求める（事前にinit_bell の実行が必要）
計算量：O(K log N)
*/
long long bell(long long N, long long K) {
long long ret = 0;
for (long long i = 0; i <= K; i++) {
ret += (pow(i, N) * fact_inv[i]) % MOD * pre_sum[K - i] % MOD;
ret %= MOD;
}
return ret;
}
int main() {
init_bell(10000000);
long long N, K;
cin >> N >> K;
cout << bell(N, K) << endl;
return 0;
}