プログラミング原論・第三章
三章は、結合演算子opを引数として受け取り、opに対するa^nを計算するアルゴリズムを最適化していきます。結合規則を意識しない左右べき乗からはじめて、共通部分式の除去・末尾再帰・ループ処理への変換・ループ式の最適化・蓄積変数の除去・特別ケース手続きへと、まあ見事に最適化していきます。コードはHaskellで書きながら読み進めているのでループ式のあたりは軽く流しましたが蓄積変数の除去は今後も使えそうなテクニックです。
power_accumulate_1 :: Num a => (a -> a -> a) -> a -> Integer -> a -> a
power_accumulate_1 op x n r
| n == 0 = r
| n == 1 = op r x
| otherwise = power_accumulate_1 op (op x x) (n `div` 2) r'
where
r' = if odd n then op r x else rこの関数をそのまま使うとすると、
*Main >power_accumulate_1 (+) 2 2 0
6
*Main >power_accumulate_1 (*) 2 2 1
8のように、opが(+)の場合と(*)の場合で初期値を使い分ける必要があります。これを
power :: Num a => (a -> a -> a) -> a -> Integer -> a
power op a n = power_accumulate_1 op a (n - 1) aのようにラップすることでより使いやすくなるわけです。
この式から、こういう実装なんて僕の脳みそではちょっと発想できそうもありません。
そして最後はこのアルゴリズムを利用してフィボナッチ数列を求めます。フィボナッチ数列は線形回帰関数で生成できるので、ベクトルの内積を使って、古い値{X(n-1),...,X(n-k+1)}を正しい位置にシフトしながら新たなXnを計算することができるということらしいのですが、この辺からはチョット理解が怪しいです。とりあえず、ベクトルの積は
で、上段にあたる(a b)や(p q)は下段(c d) (r s)から計算できるので、下段のpair<T T>を(n-k+1)乗した結果に初期値(0,1)を乗ずると答えが出るらしい(汗)です。pair<T T>(x0 x1)とpair<T T>(y0 y1)が与えられた場合、行列全体の乗算は、
となるので、この下段だけを求める演算子opは
fibonacci_matrix_multiply :: Num a => (a,a) -> (a,a) -> (a,a)
fibonacci_matrix_multiply x y = (x',y')
where
x0 = fst x
x1 = snd x
y0 = fst y
y1 = snd y
x' = x0 * (y1 + y0) + x1 * y0
y' = (x0 * y0) + (x1 * y1) となります。後は、powerにopとしてこの関数を渡せばいいんですが、Haskell作成したpower_accumulate_1はこの演算子は受け付けませんから(二項演算子はNum a => (a -> a -> a)型)タプルを受け付ける専用の関数を定義します。
power_fib_accumulate :: Num a => ((a,a)->(a,a)->(a,a))->(a,a)-> Integer->(a,a)->(a,a)
power_fib_accumulate op x n r
| n == 0 = r
| n == 1 = op r x
| otherwise = power_fib_accumulate op (op x x) (n `div` 2) r'
where
r' = if odd n then op r x else r
power_fib :: Num a => ((a,a) -> (a,a) -> (a,a))->(a,a)->Integer->(a,a)
power_fib op a n = power_fib_accumulate op a (n - 1) aで、最後に本体のfibonacciを作成して完成。
fibonacci :: Num a => Integer -> a
fibonacci n
| n == 0 = 0
| otherwise = fst r
where r = power_fib (fibonacci_matrix_multiply) (1,0) n一応、動作します。
*Main> map fibonacci [0..20]
[0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,6765]結論の節に
数学と抽象的アルゴリズムを組み合わせることで、線形回帰のn番目の要素を対する時間で生成するといった、驚くべきアルゴリズムが導き出せます。
とありますが、常にプログラムの計算量を意識してコードを書くってのはアタマではわかっているんですけどねぇー。ちょっとレベル高すぎますっって。




Comments
powerAcc のコードが推敲前のものになっていました。本当に書きたかったのは次のコードです。
powerAcc :: (Integral a) => (t -> t -> t) -> t -> a -> t -> t powerAcc _ _ 0 r = r powerAcc op x n r = powerAcc op (op x x) (div n 2) r' where r' = if odd n then op r x else rpower_accumulate_1 の型をもっと大きな範囲に適応できるようにしておけば、 fibonacci_matrix_multiply を直接渡せるようになると思います。
それから、個人的には power_accumulate_1 と fibonacci_matrix_multiply のコードは下記の powerAcc と fibMatrix の様にしたほうが分りやすい気がします。
powerAcc :: (Integral a) => (t -> t -> t) -> t -> a -> t -> t powerAcc _ _ 0 r = r powerAcc op x n r | even n = powerAcc op (op x x) (div n 2) r | otherwise = powerAcc op x (n - 1) (op r x) power' :: (Integral a) => (t -> t -> t) -> t -> a -> t power' op a n = powerAcc op a (n - 1) a fibMatrix :: Num a => (a, a) -> (a, a) -> (a, a) fibMatrix (x0, x1) (y0, y1) = (x0 * (y1 + y0) + x1 * y0, x0 * y0 + x1 * y1) fibonacci' 0 = 0 fibonacci' n = fst $ power' (fibMatrix) (1, 0) n