みなさんは一度くらい乗算・減算の使用が制限された環境でプログラミングをしなくちゃいけなかった経験、ありますよね？ないですか、そうですか。まぁそこはあるとしておきましょう。今後いつ来るか分からないですしね！
今回はそんな環境で コンビネーションの計算 を行わなければいけない状況に追い込まれた人のために、効率の良い計算方法をお伝えします。（オーバーフロー？気にしません）
ここでいうコンビネーションは、高校数学で必ず習う
{}_{n}C{}_{r} \ とか \begin{pmatrix}n \\ r\end{pmatrix}
つまり二項係数のことです。
- 実行環境
- Python 3.5
- ideone.com (2018/03/02 現在)
乗算・除算を加算・減算で再現する手法（１）
まず最初に思いつく方法ですね。
コンビネーションの式を
\begin{pmatrix}n \\ r\end{pmatrix} = \frac{n!}{r!(n-r)!}
と定義通りに表したとき、この計算に必要なのは
- $n!$ の階乗計算
- $r!$ の階乗計算
- $(n-r)!$ の階乗計算
- $r!$ と $(n-r)!$ の乗算
- $n!$ と $r!(n-r)!$ の除算
の５つの計算になります。
乗算・除算は加算・減算の組合わせ、階乗計算は乗算の組み合わせで実現できるのは周知の事実だと思いますが、いくつかの階乗計算が含まれていてなんだかたくさん計算しなくちゃいけない予感がします。
乗算、除算、階乗計算を加算・減算で愚直に実装すると
def mul(m, n):
if m > n: m, n = n, m
sum = 0
for i in range(m):
sum += n
return sum
def dev(m, n):
count = 0
while m >= n:
m -= n
count += 1
return count
def fact(n, r=1):
if n <= 1 or n <= r:
return 1
m = 1
for i in range(r, n+1):
m = mul(m, i)
return m
となり、コンビネーションの計算は
def comb1(n, r):
return dev(fact(n), mul(fact(r), fact(n-r)))
となりますね。
乗算・除算を加算・減算で再現する手法（２）
さきほどとは別の定義式として次のようなものもありますよね。
\begin{pmatrix}n \\ r\end{pmatrix} = \frac{n(n-1) \cdots (n-r+1)}{r!}
この式に必要な計算は
- $n(n-1) \cdots (n-r+1)$ の乗算
- $r!$ の階乗計算
- $n(n-1) \cdots (n-r+1)$ と $r!$ の除算
の3つとなり、さきほどより計算回数は少なくて済みそうです。
さきほどと同様に、コンビネーションの計算は次のようになりますね。
def comb2(n, r):
return dev(fact(n, n-r+1), fact(r))
動的計画法を用いた手法（逐次計算バージョン）
コンビネーションの定義式を別の視点から見てみると、もっと効率の良い計算方法があることに気づきます。
式を次のように変形してみます。
\begin{pmatrix}n \\ r\end{pmatrix} = \begin{pmatrix}n-1 \\ r\end{pmatrix} + \begin{pmatrix}n-1 \\ r-1\end{pmatrix}
ここで $\begin{pmatrix}n \ r \end{pmatrix}$ を格子点の座標と考え、$\begin{pmatrix}n \ r\end{pmatrix}$ は $\begin{pmatrix}n-1 \ r\end{pmatrix}$ と $\begin{pmatrix}n-1 \ r-1\end{pmatrix}$ に依存していると考えると、例えば $\begin{pmatrix}n \ r\end{pmatrix} = \begin{pmatrix}10 \ 5\end{pmatrix}$ は $\begin{pmatrix}9 \ 5\end{pmatrix}$ と $\begin{pmatrix}9 \ 4\end{pmatrix}$ に依存していると考えることができます。同様に依存関係を辿って行くと、次のように図示できます。もちろん依存関係の終端 $\begin{pmatrix}0 \ 0\end{pmatrix}$ は $1$ です。
あとはこの依存関係をそのまま再帰式としてプログラムを書いても良いのですが、ここは式を右から左に見る、つまり 「$\begin{pmatrix}n-1 \ r\end{pmatrix}$ と $\begin{pmatrix}n-1 \ r-1\end{pmatrix}$ が計算済みであれば $\begin{pmatrix}n \ r\end{pmatrix}$ が計算できる」と考えてみましょう。すると先の依存関係の図は次のように書き換えることができます。
さらにこの平行四辺形を長方形に整えてみましょう。
これで動的計画法でコンビネーションの計算ができる準備が整いました。プログラムは次のようになりますね。
def comb3(n, r):
r, n = r+1, n-r+1
c = np.zeros((r, n), dtype=np.int32)
c[0], c[:] = 1, 1
for i in range(1, r):
for j in range(1, n):
c[i][j] = c[i][j-1] + c[i-1][j]
return c[r-1][n-1]
これで $(r-1)(n-r-1)$ 回の加算のみ（インデックスの計算を除く）でコンビネーションの計算が可能となりました。（というかよく見るとパスカルの三角形を計算してるだけでしたね）
動的計画法を用いた方法（並列計算バージョン）
上記のやり方でも良いのですが、平行四辺形の依存関係図に立ち帰ってみると$r$軸に平行な格子点はそれぞれ同時に、つまり並列に計算できることが分かります。例えば次の図の赤で囲った部分（$n=4$ の格子点）は一つ前の格子点（$n=3$ の格子点）が全て計算済であれば並列に計算できます。
また、計算を簡略化するために、あえて $n$軸、$r$軸、直線 $n=10$、直線 $r=5$ で囲まれた部分の格子点全てを計算してみます。もちろん $\begin{pmatrix}0 \ 1\end{pmatrix}$ から $\begin{pmatrix}0 \ r\end{pmatrix}$ は全て $0$ です。
これにより計算を簡略化できるだけでなく、次のように依存関係を書き換えることができ1次元配列のみでコンビネーションの計算が行えるので、メモリの節約にもなります。
余談
上のような、単純な処理を行うノードを多数配置して相互に接続しデータを授受する機構は、シストリックアレイ（Systolic Array）とも呼ばれ、GoogleのTPU や流体力学をシミュレートするステンシル計算などに用いられています。
以上を踏まえてプログラムを書くと次のようになります（numpyの機能によりc[:r] += c[1:r+1]の部分で並列計算しているように見えますが、中ではループを回しているので実は並列計算していません）。
def comb4(n, r):
c = np.zeros(r+1, dtype=np.int32)
c[-1] = 1
for i in range(n): c[:r] += c[1:r+1]
return c[0]
ベンチマーク
「乗算・除算を加算・減算で再現する手法」2つと「動的計画法を用いた方法」2つの計4種類の実行速度を計ってみます。それぞれ8,000回実行したときの経過時間を比較してみましょう。【ベンチマークスクリプト】
| 手法 | 実行時間 [sec] |
|---|---|
| 乗算・除算を加算・減算で再現する手法（１） | 1.451 |
| 乗算・除算を加算・減算で再現する手法（２） | 0.916 |
| 動的計画法を用いた手法（逐次計算バージョン） | 1.336 |
| 動的計画法を用いた方法（並列計算バージョン） | 1.138 |
こう見ると意外と「乗算・除算を加算・減算で再現する手法（１）」と「動的計画法を用いた手法（逐次計算バージョン）」との差は開いてないですね。
一番早かったのが「乗算・除算を加算・減算で再現する手法（２）」で「動的計画法を用いた方法（並列計算バージョン）」より早いのには驚きました。計算を初期化するたびに配列を確保しているため、今回のように規模が小さいプログラムだと配列の確保がオーバーヘッドとなり計算速度が思ったよりでなかったのではないかと考えられます。
きちんと並列処理をすれば1もっと早くなると思うのですが、それは今後の課題としておきます。
最後に
今回は４種類の方法でコンビネーションの計算を行いました。他に「もっと効率の良い方法があるよ！」という方はぜひコメントお待ちしています！
-
HDLのシミュレータで厳密なクロック数を計りたいですね。 ↩