前回の記事では,OCamlという言語の変数や関数について簡単に説明した。変数を定義する構文は次の通りだった。

let 変数名 = 式

 このように,変数の「定義」は可能だが,変数への「代入」は存在しない。また,関数を定義する構文は次の通りだった。

let 関数名 引数名1 引数名2 ... 引数名n = 式

 ただし,OCamlでは,再帰関数を定義する場合は,次のようにrecを付けなければならない。

let rec 関数名 引数名1 ... 引数名n = 式

 さて,前回の記事では,べき乗を計算するpower関数の例を紹介した。power関数は,フィボナッチ数列を計算するfib関数と並んで,関数型言語のプログラムとして,最もよく引き合いに出される例の一つである。ただ,あまりにも数学的な例ばかりなので,「関数型言語は実用的ではない」という印象の原因の一つになっているのではないか,という気もする。本当は,再帰の考え方は現実の問題に対しても大いに役立つのだが,関数型言語をできるだけ簡潔に説明するという目的では,やはり数学の例題も便利だ。そこで,今回はpower関数とfib関数について,ちょっとマニアックに掘り下げてみたい。

power関数の改良

 前回は,整数mのn乗を計算する関数powerを,以下のように定義した。

let rec power m n =
  if n <= 0 then 1 else
  power m (n - 1) * m

 しかし,これではmやnが大きいと,すぐに整数オーバーフローが起こってしまう。32ビットCPUでは,OCamlのintは31ビットで表現されるためである(1ビットはガーベジ・コレクションのための「タグ」として使われる)。

# let rec power m n =
    if n <= 0 then 1 else
    power m (n - 1) * m ;;
val power : int -> int -> int = <fun>
# power 10 10 ;;
- : int = -737418240

 余談だが,整数オーバーフローはバグの原因となる場合もあるので,Standard MLでは仕様により例外を発生することが要求されている(OCamlが整数オーバーフローを検出しないのは,おそらく速度を重視したためと思われる)。また,Schemeでは,デフォルトで多倍長整数を使用する実装がほとんどで,整数オーバーフローはほとんど発生しない。

 OCamlにも多倍長整数のライブラリはあるが,ここでは説明を簡単にするために多倍長整数は利用せず,mを浮動小数にしよう(nは整数のままであることに注意してほしい)。

# let rec power m n =
    if n <= 0 then 1.0 else
    power m (n - 1) *. m ;;
val power : float -> int -> float = <fun>
# power 10.0 10 ;;
- : float = 10000000000.

 これでたしかに整数オーバーフローは起こらなくなった。しかし,このpower関数では,引数mとnが与えられると,計n回のかけ算が実行される。大げさな言い方をすると,計算量がO(n)すなわちnのオーダーになっている。これを減らすことはできないだろうか?

 一般に,アルゴリズムの計算量を減らすための考え方としては,「分割統治」が有名だ。問題を半分にわけて,それぞれを解いてから答えを組み合わせる,という考え方である。

 ここでも,nを半分にするという考え方が有効だ。具体的には,nが偶数の場合は「mのn乗」と「m*mのn/2乗」が等しいことを使って,次のようにpower関数を書くことができる。

let power m n =
  if n <= 0 then 1.0 else
  if n mod 2 = 0 then power (m *. m) (n / 2) else
  power m (n - 1) *. m

 こうすると,power関数によるかけ算の回数はlog(n)のオーダーになる。

 さらに,前回と同様に,「計算結果aも引数にしていまう」というテクニックで,このpower関数を末尾再帰化してみよう。

let rec power_loop m n a =
  if n <= 0 then a else
  if n mod 2 = 0 then power_loop (m *. m) (n / 2) a else
  power_loop m (n - 1) (m *. a)

let power m n = power_loop m n 1.0

 power_loopは「mのn乗にaをかけて返す関数」になっている。末尾再帰によりnが減っていき,nが0になったら最終結果としてaを返す。これも前回と同様だ。違うのは「if n mod 2 = 0 then ...」という行が加わった点だけである。

 なお,命令型言語をご存じの方は,これと同じアルゴリズムを再帰ではなくC言語のループで書いたらどうなるか,考えてみるとおもしろいかもしれない。解答の一例を最後に掲載しておく。

fib関数とその改良

 数学的な例をもう一つ考えてみよう。フィボナッチ数列を計算する関数fibだ。以下のように定義できる。

let rec fib n =
  if n <= 1 then n else
  fib (n - 1) + fib (n - 2)

 この定義では,与えられたnに対して,fib関数が何度も繰り返して呼び出される(指数オーダーの回数)。例えばfib 30を計算してみると,かなり時間がかかることがわかる(ocamloptコマンドでネイティブコードにコンパイルして実行すれば,もう少し速くなるが)。

# let rec fib n =
    if n <= 1 then n else
    fib (n - 1) + fib (n - 2) ;;
val fib : int -> int = <fun>
# let start = Sys.time ()    (* 開始時刻 *)
  let fib40 = fib 40
  let finish = Sys.time ()   (* 終了時刻 *)
  let time = finish -. start (* 計算時間 *) ;;
val start : float = 15.95
val fib40 : int = 102334155
val finish : float = 28.85
val time : float = 12.9000000000000021

 こんなに遅いのは,再帰の過程で同じnに対してfib nを何度も計算しているためだ。例えば,fib 5を計算するときにはfib 4とfib 3が呼び出されるが,fib 4を計算するときにもfib 3が呼び出される,という具合だ。

 これを改良するためには,末尾再帰版powerのときと同様に,引数を増やせば良い。具体的には,次のように書ける。

let rec fib_loop n a b =
  if n <= 0 then a else
  fib_loop (n - 1) b (a + b)
let fib n = fib_loop n 0 1

 言葉で説明するよりも,n = 5程度の小さな値について,fib nの評価を自分でシミュレートしてみるとわかりやすいだろう。式の値が等しいことを=で書くとすれば,

    fib 5
  = fib_loop 5 0 1
  = fib_loop 4 1 1
  = fib_loop 3 1 2
  = fib_loop 2 2 3
  = fib_loop 1 3 5
  = fib_loop 0 5 8
  = 5
という感じだ。

 あえて数学的に言えば,fib_loopは「nとfib(m)とfib(m+1)を受け取ると,fib(m+n)の値を返す」という関数になっている。これはnについての数学的帰納法で証明できる。もっとも,そんなことがわからなくても,関数型言語を使うことはできるが。

fib関数のさらなる改良

 以上のように改良されたfib関数は,n回の再帰でfib(n)の値を計算してくれる。では,これをさらに速くすることはできるだろうか?

 ここで,高校数学の数列の問題を思い出していただきたい。fib(n) = fib(n-1) + fib(n-2)という3項間漸化式を解くときに,ベクトルF(n)と行列Aを以下のようにおくと,F(n) = A F(n-1)と表せる,という話を聞いたことはないだろうか。

  F(n)  =  (fib(n+1), fib(n))    (2次元の縦ベクトルだが表記の都合で横に書く)
  A  =  ((1, 1), (1, 0))         (第1行が(1, 1)で第2行が(1, 0)の行列)

 すると,F(n) = A F(n-1) = A A F(n-2) = A A A F(n-3) = ...となるから,F(n)を計算するためには,Aのn乗を計算して,それに右からF(0)をかければ良いことになる。

 ここで,Aのn乗は,先のpowerと同様の方法により,オーダーlog(n)の時間で計算できる。したがって,fib nも,オーダーlog(n)の時間で計算できることになる。

 以上のアイディアをOCamlで実装すると,以下のようなプログラムが完成する(ただし,powerのときと同じく,整数オーバーフローを回避するために浮動小数を使っている)。

let rec fib_loop x y z n a b c =
  if n <= 0 then b else
  if n mod 2 = 0 then
    fib_loop
      (x *. x +. y *. y)
      (x *. y +. y *. z)
      (y *. y +. z *. z)
      (n / 2)
      a b c
  else
    fib_loop
      x y z
      (n - 1)
      (a *. x +. b *. y)
      (b *. x +. c *. y)
      (b *. y +. c *. z)

let fib n =
  fib_loop
    1.0 1.0 0.0
    n
    1.0 0.0 1.0

 引数が多いので複雑に見えるが,関数fib_loopは,((x, y), (y, z))という行列をn乗し,((a, b), (b, c))という行列に左からかけて,その2行1列の要素を返す関数である。先のpower_loopでいうところのmがx, y, zに,aがa, b, c になっただけで,関数の構造は全く同じだ。文章で説明するとややこしいが,興味のある方は紙と鉛筆で行列の計算をしてみればよくわかると思う。

 このfibならば,nが40だろうと1000だろうと,アッと言う間に計算が終了する。

# fib 40 ;;
- : float = 102334155.
# fib 1000 ;;
- : float = 4.34665576869374366e+208

 まだ余裕のある方は,このfibを再帰ではなくC言語などのループで書いたらどうなるか,考えてみていただきたい(解答の一例を最後に掲載する)。また,計算量のオーダーは変わらないが,x*y+y*zをy*(x+z)にするなど,まだ細かい最適化ができるところもある。

 今回はプログラミング言語というより数学の小難しい話になってしまった。次回からは,またプログラミング言語としてのOCamlの特徴を,できるだけわかりやすく紹介していきたい。


計算量オーダーlog(n)のpower関数をC言語のループで書いた例:

double power(double m, int n) {
  double a = 1;
  while (n > 0) {
    if (n % 2) {
      a *= m;
      n--;
    } else {
      m *= m;
      n /= 2;
    }
  }
  return a;
}

計算量オーダーlog(n)のfib関数をC言語のループで書いた例:

double fib(int n) {
  double x = 1, y = 1, z = 0, a = 1, b = 0, c = 1, old_b, old_y;
  while (n > 0) {
    if (n % 2) {
      old_b = b; /* 一時変数へ退避が必要 */
      b = b * x + c * y;
      a = a * x + old_b * y;
      c = old_b * y + c * z;
      n--;
    } else {
      old_y = y; /* 同上 */
      y = x * y + y * z;
      x = x * x + old_y * old_y;
      z = old_y * old_y + z * z;
      n = n / 2;
    }
  }
  return b;
}


【2007年12月3日 ITpro注】

元のタイトルは「計算量の工夫でプログラムは劇的に速くなる」でしたが,筆者のご要望により現行のタイトルに変更いたしました。また,最後の段落も元は

 今回は,プログラミングの説明というよりは,数学の小難しい話に感じたかもしれない。ただ,時間をかけて考えるとそれほど難しい話ではない。計算量の考え方はプログラミングを行ううえでの基本なので,ぜひ身に付けていただきたいと思う。次回からはまた,プログラミング言語としてのOCamlの特徴をできるだけわかりやすく紹介していきたい。

となっていましたが,同じく筆者のご要望により現行のように変更いたしました。ご了承ください。

著者紹介 住井 英二郎

 東北大学 大学院 情報科学研究科 助教授。ICFP(関数型プログラミングについての国際会議)が主催したプログラミング・コンテストで,東京大学の大学院生とOBのチームが2位に入賞しました(関連記事)。1位と3位は米Google社員らのチームだったそうです。昨年と一昨年は,京都大学のチームが3位にランクインしています。過去には筆者が参加した東京大学のチームが優勝したこともありました。ITにおいて「日本の学生は,米国や中国やインドの学生に負けている」などと言われますが,「そんなことはない」という自信を持ってよい結果だと思います。日本の大学や企業も,こうしたコンテストに参加する学生や社員をサポートしていただければありがたいと思うのですが,いろいろと難しいようで辛いところです。