« フィボナッチ数を一般項から求めてみる。 | トップページ | お久しぶりです。 »

2013年11月19日 (火)

フィボナッチ数・いろいろ

《 はじめに … フィボナッチ数の定義 》

またまたフィボナッチ数の話題です(笑)

前回の記事にも書きましたが、フィボナッチ数の定義は次のようになります。

\begin{eqnarray*} \begin{cases} F_{0} &=& 0 \\ F_{1} &=& 1 \\ F_{n+2} &=& F_{n+1} + F_{n} \end{cases} \end{eqnarray*}

フィボナッチ数の求め方は、これまでもいろいろと考えられてきたようです。 今回は、私の知っているいくつかの方法を紹介するとともに、実行速度についても比較してみようと思います。 (速度の比較は "GHCi" 上で行なっているので、コンパイルして実行した場合とは少し違っているかもしれません)

《 再帰的プロセス 》

まずは、フィボナッチ数の定義そのままのもの。

fib1 :: Int -> Integer
fib1 0 = 0
fib1 1 = 1
fib1 n = fib1 (n - 1) + fib1 (n - 2)
*Main> mod (fib1 10) 3
1
(0.00 secs, 0 bytes)
*Main> mod (fib1 20) 3
0
(0.05 secs, 4263488 bytes)
*Main> mod (fib1 30) 3
2
(1.87 secs, 223412808 bytes)
*Main> mod (fib1 35) 3
2
(20.30 secs, 2456928716 bytes)

この関数はとても時間がかかります。\(F_{35}\) を求めるのに 20 秒以上かかってしまいました。

SICP で詳しく解説されいるのですが、関数 fib1 は呼ばれるたびに自分を二回呼び出すので、とても非効率的です。

《 反復的プロセス 》

フィボナッチ数は反復的プロセスでも求めることができます。

fib2 :: Int -> Integer
fib2 n = loop n 0 1 
  where
    loop 0 a _ = a
    loop n a b = loop (n - 1) b (a + b)

fib2 の内部関数である "loop 関数" は初期値として \(F_{0}\) と \(F_{1}\) を持ち、それをもとに次々に \(F_{2}\), \(F_{3}\) … を計算していきます。 求めるフィボナッチ数を \(F_{n}\) とすると、この "loop" は n 回繰り返されます。ですから計算量は Θ(n) ということになると思います。

*Main> mod (fib2 1000) 3
0
(0.02 secs, 2136108 bytes)
*Main> mod (fib2 10000) 3
0
(0.05 secs, 8659128 bytes)
*Main> mod (fib2 100000) 3
0
(1.58 secs, 611890040 bytes)
*Main> mod (fib2 1000000) 3
0
(336.62 secs, 45991518544 bytes)

「再帰的プロセス」に比べれば高速ですが、まだまだ遅いですね。

《 行列を使ったもの 》

Wikipedia にも載っているのですが、フィボナッチ数の漸化式は行列でも表現できます。

\[ \left( \begin{array}{ll} F_{n + 1} & F_{n} \\ F_{n} & F_{n - 1} \end{array}\right) = \left( \begin{array}{ll} 1 & 1 \\ 1 & 0 \end{array} \right) ^ {n} \]

前回の記事で使用した "power 関数" を使用すれば行列の計算方法を "power 関数" に教えてやるだけで OK です。

行列 \( \left( \begin{array}{ll} a & b \\ c & d \end{array} \right) \) を、 Haskell のタプルを使って "(a b c d)" と表現すると、行列の積を計算する「関数 f」は次のように定義できます。

f (a1, b1, c1, d1) (a2, b2, c2, d2) = (a', b', c', d')
  where
    a' = a1 * a2 + b1 * c2
    b' = a1 * b2 + b1 * d2
    c' = c1 * a2 + d1 * c2
    d' = c1 * b2 + d1 * d2

これを使ってフィボナッチ数を求める関数を定義してみます。

power :: Integral a => (t -> t -> t) -> t -> a -> t -> t
power _ _ 0 e = e
power f a n e = power f (f a a) (div n 2) (if odd n then f a e else e)

fib3 :: Int -> Integer
fib3 n = (power f [1,1,1,0] n [1,0,0,1]) !! 1

f [a1, b1, c1, d1] [a2, b2, c2, d2] = [a', b', c', d']
  where
    a' = a1 * a2 + b1 * c2
    b' = a1 * b2 + b1 * d2
    c' = c1 * a2 + d1 * c2
    d' = c1 * b2 + d1 * d2
*Main> mod (fib3 100000) 3
0
(0.03 secs, 1940120 bytes)
*Main> mod (fib3 1000000) 3
0
(0.09 secs, 7170248 bytes)
*Main> mod (fib3 10000000) 3
0
(1.09 secs, 87110184 bytes)
*Main> mod (fib3 100000000) 3
0
(14.41 secs, 917842092 bytes)

"power 関数" は「n 乗」を Θ(log n) で計算してくれるので、これまでの方法と比べてかなり高速です。

《 SICP で紹介されているアルゴリズム 》

こちらのブログ に次のような、 SICP に載っているアルゴリズムを Haskell に移植したものが紹介されていま す(関数名は変えさせてもらってます)。

fib4 :: Int -> Integer
fib4 = iter 1 0 0 1
  where
    iter a b p q count
      | count == 0 = b
      | even count = iter a b (p*p+q*q) (2*p*q+q*q) (count `div` 2)
      | otherwise  = iter (b*q+a*q+a*p) (b*p+a*q) p q (count - 1)

SICP によると、この関数は Θ(log n) で \(F_{n}\) を求めることができるそうです。

*Main> mod (fib4 100000) 3
0
(0.00 secs, 0 bytes)
*Main> mod (fib4 1000000) 3
0
(0.09 secs, 5390084 bytes)
*Main> mod (fib4 10000000) 3
0
(0.78 secs, 68803364 bytes)
*Main> mod (fib4 100000000) 3
0
(11.31 secs, 766052000 bytes)

計算量が Θ(log n) なので、速度は fib3 とほぼ同様になります。

《 Gosper and Salamin のアイデア 》

ネットでフィボナッチ数についていろいろ調べている時に、『フィボナッチ数とリュカ数の計算法』という文献を見つけました。その中に「Gosper and Salamin のアイデア」というものが載っていました。

詳しくはこの文献を読んでいただくとして、要約すると、「順序対 (a, b) と (c, d) に対して乗法を (a, b)(c, d) = (ac + ad + bc, ac + bd) と定義すると、この乗法は結合律を満し、指数法則が成り立つ。そして f ≡ (1, 0) とすると \((F_{n}, F_{n + 1})\) = f n で表わすことができる ということのようです。これは SICP に載っている fib4 のと同じ方法のように見えます。

これを自作の "power 関数" を使って定義すると、次のようになりました。

power :: Integral a => (t -> t -> t) -> t -> a -> t -> t
power _ _ 0 e = e
power f a n e = power f (f a a) (div n 2) (if odd n then f a e else e)

fib5 :: Int -> Integer
fib5 n = fst $ power gsMulti (1, 0) n (0, 1)
  where gsMulti (a, b) (c, d) = (a * (c + d) + b * c, a * c + b * d)
*main> mod (fib5 100000) 3
0
*Main> :set +s
*Main> mod (fib5 100000) 3
0
(0.02 secs, 1686056 bytes)
*Main> mod (fib5 1000000) 3
0
(0.09 secs, 6295976 bytes)
*Main> mod (fib5 10000000) 3
0
(0.81 secs, 68134876 bytes)
*Main> mod (fib5 100000000) 3
0
(10.89 secs, 727571692 bytes)

"power 関数" を使っているので、計算量はやはり Θ(log n) になります。

《 一般項の式を使ったもの 》

最後に前回の記事に書いた「一般項の式を使ったもの」です。

import Data.Ratio

power :: Integral a => (t -> t -> t) -> t -> a -> t -> t
power _ _ 0 e = e
power f a n e = power f (f a a) (div n 2) (if odd n then f a e else e)

fib6 :: Int -> Integer
fib6 n = truncate . (* 2) . snd $ power f (1 % 2, 1 % 2) n (1, 0)
  where f (a1, b1) (a2, b2) = (a1 * a2 + 5 * b1 * b2, a1 * b2 + a2 * b1)
*Main> mod (fib6 100000) 3
0
(0.03 secs, 2389868 bytes)
*Main> mod (fib6 1000000) 3
0
(0.11 secs, 8976944 bytes)
*Main> mod (fib6 10000000) 3
0
(0.94 secs, 106752416 bytes)
*Main> mod (fib6 100000000) 3
0
(12.31 secs, 1174476528 bytes)

"power 関数" を使っているので、計算量は Θ(log n) のはずですが、有理数の計算をしてる分、"fib4" や "fib5" よりは少し遅くなるようです。

《 まとめ 》

これまでの結果をまとめてみます。

関数n時間 (sec)
fib13520.30
fib21,000,000336.62
fib3100,000,00014.41
fib4100,000,00011.31
fib5100,000,00010.89
fib6100,000,00012.31

今回はメモ化などの手法は考えなかったので、結果は「計算量の少なさ=計算速度」ということになっています。 fib3, fib4, fib5, fib6 はいずれも計算量が Θ(log n) なので、かなり速くフィボナッチ数を求めるこができます。

« フィボナッチ数を一般項から求めてみる。 | トップページ | お久しぶりです。 »

Haskell」カテゴリの記事

数学」カテゴリの記事

コメント

コメントを書く

コメントは記事投稿者が公開するまで表示されません。

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/112020/58605284

この記事へのトラックバック一覧です: フィボナッチ数・いろいろ:

« フィボナッチ数を一般項から求めてみる。 | トップページ | お久しぶりです。 »

2016年7月
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            
フォト

最近のトラックバック

無料ブログはココログ