掛け算と平方乗算を実装する

今日の料理

概要

『その数式、プログラミングできますか？』を読んでいたところ、一番最初に掛け算を実装する話が出てきた。掛け算は確かに身近な計算であるけれども、ではこれを効率良く計算する方法とはいったいどういうものなのかを考えるのは面白く、また再帰を書く練習になることに 気がついたので、これを実装しようと思う。今回は、言語としてOCamlを利用する。

はじめに

『その数式、プログラミングできますか？』をまったりと読んでいだところ、最初のアルゴリズムとして、掛け算の実装方法が乗っていた。掛け算はあまりにも自明なため、あまり意識しないところだが、これを再実装するのは面白い話題であることに気がついた。

まず最初に、掛け算というものを定義するとするならば、直感的には「足し算をある数で数回行なうもの」と言うことができる。もし、これを再帰的な関数として定義するとするならば、以下のようになるだろう:

let rec mul1 x n =
  if n == 0 then 0
  else if n == 1 then x
  else x + (mul1 x (n - 1))

とりあえず、話を簡単にするために、この時には負の整数の掛け算は考慮に入れないことにしよう。この方法は、例えば、掛け算がついていない電卓で、プラスのボタンを押す方法に近い。同じ数のボタンを何度も叩く、というようにだ。

しかし、これは直感的ではあるものの、愚直な方法でもある。上記の本によれば、もっと効率的な掛け算の方法として、「ロシア農民のアルゴリズム」という方法が紹介されている。これを書き下すと、下のようなコードになる。

let is_odd x = (x land 1 == 1)
let rec mul_lusia x n =
  if n == 0 then 0
  else n == 1 then x
  else
    (mul_lusia x (n lsr 1))
    + (mul_lusia x (n lsr 1))
    + (if is_odd n then x else 0)

これは、といった、結合法則をうまく利用したものである。自分は、最初、この法則を聞いてあまりにも自明すぎるため、何のためなのかさっぱりわからなかったけれども、このように、ある計算の順序を入れかえることによって、効率良く計算することができる。

余談: 奇数判定と数を半分にする方法

余談だけれども、元の本においては、ビットを利用した判定方法を利用している。それはビット論理積と、ビット右シフトという方法だ。

通常、機械においては、数というのはビットであらわされる。ビットとは、普段良く使われる十進法を二進法で表す数(binary digitの略)と言える。例えば、9の場合、1001となる。

さて、まず奇数判定のほうだが、普通、桁というのは位が上がるにつれて、その乗算となる。十進法の場合、三桁目はだし、五桁目はといったようになる。

奇数の定義を思い出すと、と書ける。二進法の場合、二桁以上の表記は全て2nとなる。従って、一桁目にビットが立っているか否かを表現すればいいことになる。ビット論理積は、ある数のビット表記を比べて、どちらもビットが立っていればいいということになる。従って、

let is_odd x = (x land 1 == 1)

と書くことができる。

今度は半分にする方法だが、先程のビットから考えるならば、桁が下がるにつれて、その桁も半分になることから、単純に右シフトすればいいということになる。ただし、考察してみればわかるように、最後の桁の1かどうかは無視されることになる。

掛け算の末尾再帰

さて、もう少し押し進めて考える。単純な再帰の場合、関数のスタック数が問題になることがあるので、思いきって末尾最適化したいと考える。

let mul n a =
  let rec mul_tailcall r n a =
    if n == 1 then r + a
    else
      mul_tailcall (if is_odd n then r + a else r ) (n lsr 1) (a + a)
  in
  if a == 1 then n
  else
    mul_tailcall n (if is_odd a then a - 1 else a) n

まず、mul_tailcallのところだが、ここで気になるのは、(a + a)の部分だろう。これに関しては、次のような考察に基づくものだ。もしaが偶数であるならば、na = n/2(a + a)とした配分法則によって表現できる。しかし、数には奇数も存在する。このとき、na = (n - 1)/2(a + a) + aとすることができる。

このaの余りの部分をどう処理するか。このとき、rに対して、この余分なaをスタックするようにしていけばよい。

逐次平方に応用する

さて、これに関しては応用的に逐次平方という考え方ができる。逐次平方とは、べき乗がとできる指数法則を利用したもので、これにより、より少ないステップで計算することが可能となる。

逐次平方の再帰的な記述は古典的に良く知られていて、例えば、SICPにも覗かれる。

要点を抜きだせば、のときのnが偶数のとき、と計算し、奇数のときはと分割していくことだ。

let rec pow x n =
  if n == 0 then 1
  else if n == 1 then x
  else if n == 2 then mul x x
  else if is_odd n then mul (pow x (n - 1)) x
  else pow (pow x (n lsr 1)) 2

これは中で再帰的に関数を呼びだしているので、結果を持つようにストックしてみよう。

let pow x n =
  let rec pow_tailcall r n a =
    if n == 1 then r * a
    else if (is_odd n) then pow_tailcall r (n - 1) (mul a r)
    else pow_tailcall (mul r r) (n lsr 1) a
  in
  if n == 0 then 1
  else if n == 1 then x
  else pow_tailcall x (n - 1) x

蛇足: OUnit

ちなみに、このプログラムについでは、OUnitを使いながら実装した。テストをやってみると、意外とボロが出るんだなあと痛感した。

open OUnit2;;
open Multi;;

let test_mul1 test_ctxt =
  assert_equal 3 (mul 3 1);;
let test_mul test_ctxt =
  assert_equal 6 (mul 2 3);;
let test_pow0 test_ctxt =
  assert_equal 1 (pow 3 0);;
let test_pow1 test_ctxt =
  assert_equal 3 (pow 3 1);;
let test_pow72 test_ctxt =
  assert_equal 49 (pow 7 2);;
let test_pow75 test_ctxt =
  assert_equal 16807 (pow 7 5);;

let suite = "掛け算と逐次平方の実装" >:::
              ["1回だけかける" >:: test_mul1;
               "複数回かける" >:: test_mul;
               "xの0乗は1" >:: test_pow0;
               "xの1乗はx" >:: test_pow1;
               "7の2乗は49" >:: test_pow72;
               "7の5乗は16807" >:: test_pow75
               ];;

let () = run_test_tt_main suite;;

実行に関しては、いい方法なのか解らないが、make.shに以下のようなコマンドを書いて実行権限を与えていた。

ocamlfind ocamlc -package oUnit -linkpkg -o multi -g multi.ml multi_test.ml
./multi

OCamlのエコシステムに関しては正直よくわからないところがあるので、改めて調べなおしたほうがいいかもしれない。

おわりに

ということで、簡単な掛け算の実装をやってみたのだけれども、中学生のような掛け算の法則をふりかえることで、いろいろと考えることができたので、簡単なところからやっていくのも重要だなあと思ったりした。

その数式、プログラムできますか？