第2種スターリング数を求めるアルゴリズム
第2種スターリング数は「写像12相」と呼ばれるものを通して学ぶと、他の数え上げ問題との関わりが分かり全体像がスッキリします。ぜひ確認してみることをおすすめします。
「写像12相」で典型的な数え上げ問題のパターン総整理
第2種スターリング数とは
第2種スターリング数 \(S(n,k)\) は、以下の3つの性質を持つ数です。
- \(S(0,0)=1, ~ \forall i_{\geq 1}. S(i,0)=0\) として、\(S(n,k) = k S(n-1,k) + S(n-1, k-1)\) が成立する
- 区別できる \(n\) 個の要素を区別できない \(k\) 個のブロックに分割する方法の数
- \(f_k(x)=x(x-1)\cdots(x-k+1)\) として、\(x^n=\displaystyle\sum_{k=1}^nS(n,k)f_k(x)\) が成立する
これら3つは全て同値な性質であり、1つを定義とすれば他の2つを導出することができます。
アルゴリズム
1つ目の性質である \(S(n,k) = k S(n-1,k) + S(n-1, k-1)\) を用いることで、動的計画法によって計算することができます。
もしくはスターリング数に関する以下の有名公式を利用することもできます。
$$ S(n,k) = \frac{1}{k!}\sum_{i=0}^{k} (-1)^{k-i} {}_{k}{\mathrm C}_i i^n $$
二項係数の計算ができれば、この右辺を用いて計算することが可能です。
動的計画法によるアルゴリズム：
- 二次元配列 S を作り0 で初期化する
- S(0,0)=1 とする
- 有名公式 \(S(n,k) = k S(n-1,k) + S(n-1, k-1)\) を利用して、動的計画法により順番に \(S(n,k)\) を計算する。
※ 以下のようなテーブルを作成することになります。
| n ＼ k | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | ||||||||
| 1 | 0 | 1 | |||||||
| 2 | 0 | 1 | 1 | ||||||
| 3 | 0 | 1 | 3 | 1 | |||||
| 4 | 0 | 1 | 7 | 6 | 1 | ||||
| 5 | 0 | 1 | 15 | 25 | 10 | 1 | |||
| 6 | 0 | 1 | 31 | 90 | 65 | 15 | 1 | ||
| 7 | 0 | 1 | 63 | 301 | 350 | 140 | 21 | 1 | |
| 8 | 0 | 1 | 127 | 966 | 1701 | 1050 | 266 | 28 | 1 |
二項係数によるアルゴリズム：
- \(\frac{1}{k!}\sum_{i=0}^{k} (-1)^{k-i} {}_{k}{\mathrm C}_i i^n\) を計算する
※ 有名公式 \(S(n,k) = \frac{1}{k!}\sum_{i=0}^{k} (-1)^{k-i} {}_{k}{\mathrm C}_i i^n\) の右辺を計算しています。
※ 二項係数は前計算を行っておくことで、素早く計算することができます。詳細は以下の記事を御覧下さい。
競プロでよく使う二項係数(nCk)を素数(p)で割った余りの計算と逆元のまとめ
※ \(i^n\) に関しては、繰り返し2乗法などを用いることで高速に計算ができます。
繰り返し二乗法によるべき乗(pow(x,n))の計算のアルゴリズム
計算量
- 動的計策法によるアルゴリズム： \(O(kn)\)
- 二項係数によるアルゴリズム： \(O(k \log n)\)
動的計策法によるアルゴリズムでは、\(S(n,k)\) を求めるため動的計画法により \(n \times k\) のテーブルを更新することになります。1つのマスを計算するのに \(O(1)\) かかるので、全体では \(O(nk)\) だけかかってしまいます。
\(\frac{1}{k!}\sum_{i=0}^{k} (-1)^{k-i} {}_{k}{\mathrm C}_i i^n\) を求めるアルゴリズムでは、二項係数については前計算を行うことで \(O(1)\) 程度で計算できます。\(i^n\) を求めるのには\(O(\log n)\) だけかかるので、\(k\) 回分計算すると全体では \(O(k \log n)\) になります。
C++での実装例
AOJ DPL_5_I Balls and Boxes 9 を解く実装例です。第2種スターリング数を求めます。
答えは非常に大きくなることがあるので、 素数 1000000007 で割った数を答えにしています。
動的計画法によるアルゴリズム
\(n \times k\) のテーブルを更新します。複数の第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;
}
}
cout << S[N][K] << endl;
return 0;
}
二項係数によるアルゴリズム
二項係数の計算については、以下の記事を御覧下さい。
競プロでよく使う二項係数(nCk)を素数(p)で割った余りの計算と逆元のまとめ
べき乗の計算に関しては、繰り返し2乗法などを用いることで高速に計算ができます。
繰り返し二乗法によるべき乗(pow(x,n))の計算のアルゴリズム
#include <bits/stdc++.h>
using namespace std;
const int MOD = 1000000007; // 素数
vector<long long> fact, fact_inv, inv;
/* init_nCk :二項係数のための前処理
計算量:O(n)
*/
void init_nCk(int SIZE) {
fact.resize(SIZE + 5);
fact_inv.resize(SIZE + 5);
inv.resize(SIZE + 5);
fact[0] = fact[1] = 1;
fact_inv[0] = fact_inv[1] = 1;
inv[1] = 1;
for (int i = 2; i < SIZE + 5; i++) {
fact[i] = fact[i - 1] * i % MOD;
inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD;
fact_inv[i] = fact_inv[i - 1] * inv[i] % MOD;
}
}
/* nCk :MODでの二項係数を求める(前処理 int_nCk が必要)
計算量:O(1)
*/
long long nCk(int n, int k) {
assert(!(n < k));
assert(!(n < 0 || k < 0));
return fact[n] * (fact_inv[k] * fact_inv[n - k] % MOD) % 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;
}
long long S(long long n, long long k) {
long long ret = 0;
for (long long i = 0; i <= k; i++) {
if ((k - i) % 2 == 0) {
ret = ret + nCk(k, i) * pow(i, n) % MOD;
ret %= MOD;
} else {
ret = ret + (MOD - nCk(k, i) * pow(i, n) % MOD);
ret %= MOD;
}
}
ret = ret * fact_inv[k] % MOD;
return ret;
}
int main() {
init_nCk(10000000);
long long N, K;
cin >> N >> K;
cout << S(N, K) << endl;
return 0;
}