組み合わせ数（二項係数）を計算する

「○個のものから□個だけ取り出す」という処理はときどき必要となりますが、それを行った場合の「数」だけ考える場合も、地味に面倒だったりします。

定義とその性質

このような、「$n$ 個のものから $k$ 個だけ取り出す場合の数」は、二項定理の係数としても現れることから、二項係数という呼び名があります。いろいろな表記法がありますが、ここでは $_n\mathrm{C}_k$ と書いてみることにします。

まず、掛け算として書くと、

_n\mathrm{C}_k=\frac{n!}{k!(n-k)!} = \frac{n\times(n-1)\times\dots\times(n-k+1)}{k!}\tag{1}

のように書けます。また、この二項係数はいくつか満たす性質があります。

\begin{align}
_n\mathrm{C}_k&={}_n\mathrm{C}_{n-k}\tag{2}\\
  &= {}_{n-1}\mathrm{C}_k + {}_{n-1}\mathrm{C}_{k-1}\qquad(k \gt 0)\tag{3}\\
_n\mathrm{C}_0&=1\tag{4}\\
_n\mathrm{C}_1&=n\tag{5}\\
\end{align}

とりわけ(3)は、パスカルの三角形として知られています。

Rubyの場合…何も考えなくてもじゅうぶん速い

Rubyの場合、こんな場合に便利な言語仕様が1つあって、それは「整数演算で値が大きくなっても桁あふれせず、勝手に多倍長整数（Bignum）になってくれる」ということです。ということで、桁あふれを気にせず、いきなり(1)の式を使うことができます（ただ、項数が増えても無駄なので、$n/2<k$ の場合には (2)を適用して、掛ける項数を減らすことだけはしています）。

combi.rb

def combi(n,k)
  k=n-k if 2*k > n
  return 1 if k == 0
  ((n-k+1)..n).reduce(&:*)/((1..k).reduce(&:*))
end

なお、RubyではRangeもEnumerableなので、そのまま範囲内の値についてreduceできます（関連）。

そして、速度についても申し分なく、$_{1000}\mathrm{C}_{500}$ の計算が1ミリ秒で完了する程度に高速でした。

JavaScriptの場合…

JavaScriptの場合、標準で多倍長整数をサポートしておらず、(1)をそのまま書けば浮動小数点数での演算となって、2進法で53ビットの精度を越えてしまえば正確な値が出なくなります。ということで、正確な値を出すための戦略は2つ考えられます。

1. より小さいものの足し算から算出する
2. 途中で適宜約分しながら掛け算を進める

足し算で計算

(3)を使えば、より小さな値の計算へと変換することができます。ということで、まずはシンプルに実装してみましょう。

シンプルな実装

function combi_simple(n, k){
  if(k*2 > n) k=n-k;
  if(k===0) return 1;
  if(k===1) return n;
  return combi_simple(n-1, k) + combi_simple(n-1, k-1);
}

これで、解が$2^{53}$ 以下になるのでJavaScriptでも正確に表せる $_{100} \mathrm{C}_{13}$ を計算させようとしましたが…10秒しても応答がありませんでした。ただのしかばねのようだ。

キャッシュさせてみる

実は、この計算にはかなりの無駄があります。

\begin{align}
_{100} \mathrm{C}_{13}&={}_{99}\mathrm{C}_{13}+{}_{99}\mathrm{C}_{12}\\
&={}_{98}\mathrm{C}_{13}+{}_{98}\mathrm{C}_{12}+{}_{98}\mathrm{C}_{12}+{}_{98}\mathrm{C}_{11}
\end{align}

というように、展開していく中で同じ項が何度も登場するのですが、そのたびごとに計算してもあまり意味がありません。そこで、メモ化させてみることにします。

メモ化

var memo={};
function combi_memorized(n, k){
  //メモるまでもなく出る値はそのまま返す
  if(k*2 > n) k=n-k;
  if(k===0) return 1;
  if(k===1) return n;
  //メモを確認
  var memo_key = n + ',' + k;
  if(!memo[memo_key]){
    memo[memo_key] = combi_memorized(n-1, k) + combi_memorized(n-1, k-1);
  }
  return memo[memo_key];
}

こうすることで、一瞬で値を返すようになりました。

約分して掛け算

もう1つの戦略の「約分して掛け算」について考えてみましょう。できるだけ約分の回数を減らすために、「$2^{53}$ を超えない範囲で掛け算をしてから、最大公約数で割る」ことにしてみます。

約分

// 2^53より少し小さな値
var MAX=9e15;

function combi_reduction(n, k){
  //すぐに出る値はそのまま返す
  if(k*2 > n) k=n-k;
  if(k===0) return 1;
  if(k===1) return n;
  function gcd(x, y){
    var t;
    while(y){
      t = x % y;
      x = y;
      y = t; 
    }
    return x;
  }
  var num=1, denom=1, tmp, i;
  var nums=[], denoms=[];
  for(i=0; i<k; ++i){
    nums.push(n-i);
    denoms.push(k-i);
  }
  while(denoms.length){
    while(denoms.length){
      tmp = denom * denoms[denoms.length-1];
      if(tmp > MAX) break;
      denom = tmp;
      denoms.pop();
    }
    while(nums.length){
      tmp = num * nums[nums.length-1];
      if(tmp > MAX) break;
      num = tmp;
      nums.pop();
    }
    tmp = gcd(num, denom);
    num /= tmp;
    denom /= tmp;
  }
  for(i=0;i<nums.length;++i){
    num *= nums[i];
  }
  return num / denom;
}

この2つの速度比較をjsPerfに投稿しましたので確認できますが、メモ化のほうが速そうです。