第2種スターリング数を求めるアルゴリズム

2020年10月22日数学数え上げ,写像12相,スターリング数,第2種スターリング数

第2種スターリング数は「写像12相」と呼ばれるものを通して学ぶと、他の数え上げ問題との関わりが分かり全体像がスッキリします。ぜひ確認してみることをおすすめします。
「写像12相」で典型的な数え上げ問題のパターン総整理

第2種スターリング数とは

第2種スターリング数 \(S(n,k)\) は、以下の3つの性質を持つ数です。

  1. \(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)\) が成立する
  2. 区別できる \(n\) 個の要素を区別できない \(k\) 個のブロックに分割する方法の数
  3. \(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 $$

二項係数の計算ができれば、この右辺を用いて計算することが可能です。

動的計画法によるアルゴリズム:

  1. 二次元配列 S を作り0 で初期化する
  2. S(0,0)=1 とする
  3. 有名公式 \(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

二項係数によるアルゴリズム:

  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;
}