二項係数 ${}_n C_k = \binom{n}{k}$ の計算アルゴリズムを確認し, Rust でベンチマークします.

アルゴリズム

最初に注意ですが, 以下のコード例のうちいくつかは $k > n/2$ のとき $\binom{n}{k} = \binom{n}{n-k}$ となることを用いると計算量が削減できます. が, コードを読みやすくするためにこの処理を省略しています.

階乗を用いる場合

最も安直な二項係数の計算法として, 階乗を用いた表示
$$\binom{n}{k} = \frac{ n! }{ k! (n-k)! }$$
を利用するのがまず思い付きます. コードはこんな感じでしょうか.

pub fn binom_fact( n: i32, k: i32 ) -> i32 {
    factorial( n ) / factorial( k ) / factorial( n-k )
}

fn factorial( i: i32 ) -> i32 {
    if i <= 1 { 1 } else { i * factorial(i-1) }
}

ただ, 階乗はあっという間にオーバーフローするのでこれでは実用的ではないです. 例えば i32 なら $13!$ を計算することはできないため, 上のコードは $n \leq 12$ にしか使えません.

漸化式 (Pascal の三角形) を用いる場合

Pascal の三角形の $n$ 段目と $n-1$ 段目の要素の間には次の漸化式
$$\binom{n}{k} = \binom{n-1}{k-1} + \binom{n-1}{k}$$
が成り立ちます. これはすぐに再帰的なコードが書けます.

pub fn binom_rec( n: i32, k: i32 ) -> i32 {
    if k == 0 || k == n {
        1
    } else {
        binom_rec( n-1, k-1 ) + binom_rec( n-1, k )
    }
}

しかしながら, この方式では 1 段遡るごとに 2 個の要素を計算することになるので, $2^{n-1}$ 回の関数呼び出しが必要になります. これだと大きな $n$ に対して急速に遅くなります. ただ, この $2^{n-1}$ 回の評価には実際には同一の係数が何度も現れているため, かなり無駄が多いことになります. なので $m$ 段目 ($m = 1, 2, \cdots, n-1$) の二項係数を全部計算してしまえば良いのです. これは加算が $\frac{1}{2} n (n+1) = \mathcal{O} (n^2)$ 回含まれ, 必要なメモリは $\mathcal{O} (n)$ です.

pub fn binom_pascal( n: i32, k: i32 ) -> i32 {
    if k == 0 || k == n {
        return 1;
    }
    let k = k as usize;

    let mut p = vec![ 1 ];

    for _ in 0..n-1 {
        let mut c = vec![ 1 ];

        for x in p.windows(2) {
            c.push( x[0] + x[1] );
        }
        c.push( 1 );

        p = c;
    }

    p[k-1] + p[k]
}

ただし, Pascal の三角形は左右対称なのでそれを考慮すると計算量を約半分にできます (コード省略). この方法は答えがオーバーフローしない範囲に収まるならば必ずオーバーフローせずに答えを求めることができるという利点もあります.

Knuth の方法

このブログ記事によると, 二項係数の階乗による表式を
$$\binom{n}{k} = \frac{n}{1} \cdot \frac{n-1}{2} \cdot \frac{n-2}{3} \cdots \frac{n-k+1}{k}$$
と書き換える賢いアルゴリズムが Knuth の教科書に載っているそうです. これは左から順に計算していくことができます: 整数 $d$ で割るとき, 分子は $d$ 個の隣接する整数の積ですからそのどれかに必ず因子 $d$ が含まれます.

pub fn binom_knuth(n: i32, k: i32) -> i32 {
    (0..n + 1)
        .rev()
        .zip(1..k + 1)
        .fold(1, |mut r, (n, d)| { r *= n; r /= d; r })
}

これは $\mathcal{O} ( k )$ 回の積および商演算からなります. 関数の再帰呼び出しの形に書くこともできます.

同様の考え方で, 二項係数を
$$\binom{n}{k} = \frac{n-k+1}{1} \cdot \frac{n-k+2}{2} \cdot \frac{n-k+3}{3} \cdots \frac{ n-1 }{ k-1 } \cdot \frac{n}{k}$$
と書き直し, 再帰形にまとめたものがこちらのブログ記事に書かれている方法です.

pub fn binom_asc( n: i32, k: i32 ) -> i32 {
    if k == 0 || k == n {
        1
    } else {
        binom_asc( n-1, k-1 ) * n / k
    }
}

ただし答えがオーバーフローしない範囲であっても大きな $n$ の場合計算途中でオーバーフローが起きます. 上のどちらの形でも, i32 だと $\binom{n}{n/2}$ は $n = 30$ で, $\binom{n}{n/3}$ は $n = 34$ でオーバーフローしました.

ベンチマーク

それでは $\binom{29}{14} = 77,558,760$ を計算する速度を競ってもらいましょう. なお binom_pascal では対称性を利用して計算量を削減したコードを使用しています.

#![feature(test)]
extern crate test;
use test::Bencher;

const N: i32 = 29;
const K: i32 = 14;
const C: i32 = 77558760;

#[bench]
fn knuth(b: &mut Bencher) {
    b.iter(|| {
        let n = test::black_box(N);
        let k = test::black_box(K);

        assert_eq!( binom_knuth(n, k), C );
    });
}

#[bench]
fn asc(b: &mut Bencher) {
    b.iter(|| {
        let n = test::black_box(N);
        let k = test::black_box(K);

        assert_eq!( binom_asc(n, k), C );
    });
}

#[bench]
fn rec(b: &mut Bencher) {
    b.iter(|| {
        let n = test::black_box(N);
        let k = test::black_box(K);

        assert_eq!( binom_rec(n, k), C );
    });
}

#[bench]
fn pascal(b: &mut Bencher) {
    b.iter(|| {
        let n = test::black_box(N);
        let k = test::black_box(K);

        assert_eq!( binom_pascal(n, k), C );
    });
}

結果は binom_asc が 125 ns, binom_knuth が 75 ns, binom_pascal が 2,176 ns, binom_rec が 180 ms ということで, Knuth の方法が最適という結論でした. ただオーバーフローし得ることに注意が必要です.

参考文献

こちらの記事などは ${}_n C_k \ \mathrm{mod} \, p$ を計算しています.

[Rust] 二項係数の計算