べき乗計算を対数時間で行う
考えかた
たとえば3の10乗を考えてみよう。
3 * 3 * 3 * 3 * 3 * 3 * 3 * 3 * 3 * 3
だ。9回の乗算が行われている。乗算の回数を減らすために以下のようにすることができる。
9 = 3 * 3
9 * 9 * 9 * 9 * 9
こうすると5回の乗算で計算できる。同じようにさらに以下のようにする。
9 = 3 * 3
81 = 9 * 9
81 * 81 * 9
これで乗算は4回になった。aのn乗は以下のようにすると効率的に計算できる。
aのn乗を計算するのに蓄積変数としてrを用意して
- nが偶数ならばaをaの2乗にしてnを半分にする。rはそのまま
- nが奇数ならばaはそのまま。nを1減らしてrにaをかける
このときa ^ n * rはいつも同じ値となる。r = 1で始めてn = 0になるまで計算を続けるとr = a ^ nとなる。
Haskellでは
Haskellで表現すると単純なかけ算では
expt :: Integer -> Integer -> Integer
expt a n = exptIter a n 1
exptIter :: Integer -> Integer -> Integer -> Integer
exptIter _ n r | n < 1 = r
exptIter a n r = exptIter a (n - 1) (r * a)
となる。これは掛け算をn回くりかえすのでO(n)となる。かけ算の回数を減らすやりかたでは
exptIter :: Integer -> Integer -> Integer -> Integer
exptIter _ n r | n < 1 = r
exptIter a n r
| even n = exptIter (a ^ 2) (n `div` 2) r
| otherwise = exptIter a (n - 1) (r * a)
のようになる。このようにすると時間効率はO(log n)となる。変数rについてポイントフリースタイルにすると
exptIter :: Integer -> Integer -> Integer -> Integer
exptIter _ n | n < 1 = id
exptIter a n
| even n = exptIter (a ^ 2) (n `div` 2)
| otherwise = exptIter a (n - 1) . (* a)
のようになる。O(n)からO(log n)へと変換できるのは(* a)を2回行うことと(* a ^ 2)を1回行うこととが同じ結果となるからだ。
対数時間でフィボナッチ数列の第n項を求める
同様の変換をフィボナッチ数列の第n項目を求めるのにも使える。
指数時間で
まずフィボナッチ数列の第n項を定義からそのままHaskellにすると
fib :: Integer -> Integer
fib n | n < 2 = n
fib n = fib (n - 1) + fib (n - 2)
となる。これはfib nを求めるのにはfib (n - 1)とfib (n - 2)を計算し、fib (n - 1)を求めるのにはfib (n - 2)とfib (n - 3)を計算し、fib (n - 2)を求めるのにはfib (n - 3)とfib (n - 4)を計算し、といったように再帰をたどるごとに処理は2倍2倍で増えていく。
よって時間効率はO(2 ^ n)となる。
線形時間で
ここで隣りあう2項をペアとして考えることにする。
0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89 ...
を以下のようにペアにする。
(0, 1), (1, 1), (1, 2), (2, 3), (3, 5), (5, 8) ...
このようにすると次の項を求めるのに必要なのは直前の1ペアだけとなる。
fib :: Integer -> Integer
fib n = fst $ fib n (0, 1)
fibIter :: Integer -> (Integer, Integer) -> (Integer, Integer)
fibIter n r | n < 1 = r
fibIter n (a, b) = fibIter (n - 1) (b, a + b)
これで時間効率をO(n)とすることができる。
対数時間で
さて、この関数を効率的ではないべき乗の関数と比較してみよう。nを1減らして蓄積変数の値を「次の値」に変換しているのは共通している。
乗算では
n -> n - 1, r -> r * a
であり、フィボナッチ数列の場合は
n -> n - 1, (a, b) -> (b, a + b)
である。
べき乗関数でnが偶数のときについて変換のしかたを変えた。
- nが偶数のとき、a -> a ^ 2, n -> n `div` 2, r -> r
- nが奇数のとき、a -> a, n -> n - 1, r -> r * a
つまり、aを2回かける代わりにa ^ 2をかけることで乗算の回数を減らした。
フィボナッチ数列についても関数を2回適用する代わりとなる別の関数が存在すれば対数時間にすることができるはずだ。
ここでTpqという変換を考える。この変換は以下のような変換だ。
(a, b) -> (pa + qb, qa + (p + q)b)
たとえばT01はp = 0, q = 1なので
(a, b) -> (0a + 1b, 1a + (0 + 1)b)
つまり(a, b) -> (b, a + b)となる。つまりフィボナッチ数列の次の項を求める変換はTpqのうちp = 0, q = 1としたものである。
ここでTpqを2回適用したT^2pqを考えてみよう。
(a, b) -> (pa + qb, qa + (p + q)b)
とし、さらに同じ変換をするので
(pa + qb, qa + (p + q)b) ->
(p(pa + qb) + q(qa + (p + q)b),
q(pa + qb) + (p + q)(qa + (p + q)b))
となる。これを整理して
((p^2 + q^2)a + (2pq + q^2)b,
(2pq + q^2)a + ((p^2 + q^2) + (2pq + q^2))b
となる。つまり
T^2pq = T(p^2 + q^2)(2pq + q^2)
だ。これでフィボナッチ数列における整数のペアの変換を「二乗二乗」にすることができるようになった。
まずはO(n)の関数をTpqを使ったやりかたで書いてみよう。
tpq :: (Integer, Integer) -> (Integer, Integer) -> (Integer, Integer)
tpq (p, q) (a, b) = (p * a + q * b, q * a + (p + q) * b)
fib :: Integer -> Integer
fib n = fst $ tpqn (0, 1) n (0, 1)
tpqn :: (Integer, Integer) -> Integer -> (Integer, Integer) -> (Integer, Integer)
tpqn _ n r | n < 1 = r
tpqn pq n ab = tpqn pq (n - 1) (tpq pq ab)
関数tpqnのnが偶数のときにpqを「二乗」するようにする。
tpqn :: (Integer, Integer) -> Integer -> (Integer, Integer) -> (Integer, Integer)
tpqn _ n r | n < 1 = r
tpqn pq@(p, q) n ab
| even n = tpqn (p ^ 2 + q ^ 2, 2 * p * q + q ^ 2) (n `div` 2) ab
| otherwise = tpqn pq (n - 1) (tpq pq ab)
これでO(log n)のフィボナッチ関数が作成できた。関数tpqnをabに関してポイントフリースタイルにしておく。
tpqn :: (Integer, Integer) -> Integer -> (Integer, Integer) -> (Integer, Integer)
tpqn _ n | n < 1 = id
tpqn pq@(p, q) n
| even n = tpqn (p ^ 2 + q ^ 2, 2 * p * q + q ^ 2) (n `div` 2)
| otherwise = tpqn pq (n - 1) . tpq pq
Tpqをn回適用する関数を求めるのに
- nが偶数ならばT^2pqをn / 2回適用し
- nが奇数ならばTpqをn - 1回適用する関数に関数tpqを1回追加で行う
まとめ
フィボナッチ数列の次項を求める変換をTpqというより一般的な変換のp = 0, q = 1としたものと考えた。Tpqを2回適用したT^2pqがp, qをp^2 + q^2, 2pq + q^2としたT(p^2 + q^2)(2pq + q^2)とは同じである。フィボナッチ数列の変換を「2乗2乗」にしていくことができ、これによって対数時間で第n項を求めることが可能となった。
おまけ
反復的プロセスで説明してきた。再帰的プロセスでの解は以下のようになる。
tpq :: (Integer, Integer) -> (Integer, Integer) -> (Integer, Integer)
tpq (p, q) (a, b) = (p * a + q * b, q * a + (p + q) * b)
fastFib :: Integer -> Integer
fastFib = fst . fastFibRec (0, 1)
fastFibRec :: (Integer, Integer) -> Integer -> (Integer, Integer)
fastFibRec _ n | n < 1 = (0, 1)
fastFibRec pq@(p, q) n
| even n = fastFibRec (p ^ 2 + q ^ 2, 2 * p * q + q ^ 2) (n `div` 2)
| otherwise = tpq pq $ fastFibRec qp (n - 1)