大きな数字でもコンビネーションの計算を速く求めるには？

Question

以下の工夫なしのコードでは、コンビネーションの計算をするのに
１分以上かかっています。

def ncr(n, r)
  r = n - r if r > n - r
  return 0 if r < 0
  return 1 if r == 0
  return (n - r + 1..n).inject(:*) / (1..r).inject(:*)
end

N = 1000
cnt = 0
(0..N).each{|a|
  (0..a).each{|b|
    cnt += ncr(a, b)
  }
}
p cnt == 2 ** (N + 1) - 1

大きな数字でもコンビネーションの計算を
速く求められるようにするには
ncr をどのように定めればよろしいのでしょうか？

（追記）
ちなみに、Python版（工夫なし）では同じ計算を１分以内で行います。

from math import factorial

def product(iterable):
    prod = 1
    for n in iterable:
        prod *= n
    return prod

def npr(n, r):
    return product(range(n - r + 1, n + 1))

def ncr(n, r):
    if r > n // 2:
        r = n - r
    return npr(n, r) // factorial(r)

N = 1000
cnt = 0
for a in xrange(N + 1):
    for b in xrange(a + 1):
        cnt += ncr(a, b)

print cnt == 2 ** (N + 1) - 1

Answer 1

入山徳夫氏によるnCrを高速に求めるアルゴリズム(追記：リンク先変更)
人間が手計算でnCrを計算するとき、約分して分母を消して、残った分子を掛けますが、
それをそのまんま再現したようなアルゴリズムです。

もしくはパスカルの三角形を使って加算だけで求める手もあります。

どっちが速いかは知らないので早かった方をお使い下さい

あと実装は自分でやって下さい。
リンク先をRubyにただ翻訳すればいいと思います。

Comments

[manman]

入山徳夫氏によるnCrを高速に求めるアルゴリズム を試しましたが、nが4以上で誤った値をかえします。

[ozwk]

リンク先の実装が間違えてますね。申し訳ありません。
リンク先差し替えてついでに検算もしておきました。

Answer 2

入山徳夫氏によるnCrを高速に求めるアルゴリズムの
Ruby版およびPython版です。
（Ruby版およびPython版どちらの場合も、）
質問の例だと残念ながら工夫なしの方が速いのですが、
N が大きいときの計算をさせると工夫なしより速くなりました。

def ncr(n, r)
  r = [r, n - r].min
  return 1 if r == 0
  return n if r == 1
  numerator = (n - r + 1..n).to_a
  denominator = (1..r).to_a
  (2..r).each{|p|
    pivot = denominator[p - 1]
    if pivot > 1
      offset = (n - r) % p
      (p - 1).step(r - 1, p){|k|
        numerator[k - offset] /= pivot
        denominator[k] /= pivot
      }
    end
  }
  result = 1
  (0..r - 1).each{|k|
    result *= numerator[k] if numerator[k] > 1
  }
  return result
end

N = 10000
cnt = 0
(0..N).each{|a|
  cnt += ncr(N, a)
}
p cnt == 2 ** N

def ncr(n, r):
    r = min(r, n - r)
    if r == 0: return 1;
    if r == 1: return n;
    numerator = [n - r + i + 1 for i in range(r)]
    denominator  = [i + 1 for i in range(r)]
    for p in range(2, r + 1):
        pivot = denominator[p - 1]
        if pivot > 1:
            offset = (n - r) % p
            for k in range(p - 1, r, p):
                numerator[k - offset] /= pivot
                denominator[k] /= pivot
    result = 1
    for k in range(r):
        if numerator[k] > 1: 
            result *= numerator[k]
    return result

N = 10000
cnt = 0
for a in xrange(N + 1):
    cnt += ncr(N, a)

print cnt == 2 ** N