« 2012年10月 | トップページ | 2016年7月 »

2013年11月

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) なので、かなり速くフィボナッチ数を求めるこができます。

2013年11月14日 (木)

フィボナッチ数を一般項から求めてみる。

フィボナッチ数を一般項から求めてみる

《 はじめに 》

フィボナッチ数とは、

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

という式で定義される数ですが、次のような一般項を表わす式もわかっています。

\[ F_{n} = \frac{1}{\sqrt{5}} \left\{ \left( \frac{1 + \sqrt{5}}{2} \right) ^ n - \left( \frac{1 - \sqrt{5}}{2} \right) ^ n \right\} \]

今回はこの一般項の式を用いてフィボナッチ数を求める方法を考えてみました。


《 一般項のとおりに計算してみると… 》

「一般項の式があるのなら、それをそのまま使えばいいのでは…」と考える方もいると思います。
では、試してみましょう。

phi  = (1 + sqrt 5) / 2
phi' = (1 - sqrt 5) / 2

fib1 :: Int -> Integer
fib1 n = truncate $ (phi ^ n - phi' ^ n) / (sqrt 5)

この関数を "GHCi" に読み込んで実行してみます。

*Main> map fib1 [1..10]
[1,1,2,3,5,8,13,21,34,54]

F10 (10番目のフィボナッチ数) は "55" のはずなのに、"fib1 10" の結果は "54" になってしまいました。
これは fib1 が浮動小数点数で計算した値を整数に変換しているため、計算誤差が生じていると考えられます。

正しいフィボナッチ数を求めるためには、他の方法を考えないといけないようです。


《 一般項の式をよく見てみると… 》

一般項の式を使って、実際にいくつかのフィボナッチ数を計算してみましょう。

\begin{eqnarray*} & F_{1} = & \frac{1}{\sqrt{5}} \left\{ \left( \frac{1}{2} + \frac{1}{2} \sqrt{5} \right) - \left( \frac{1}{2} - \frac{1}{2} \sqrt{5} \right) \right\} & = & \frac{1}{\sqrt{5}} ・ \frac{2}{2} \sqrt{5} & = & 1 \\ & F_{2} = & \frac{1}{\sqrt{5}} \left\{ \left( \frac{6}{4} + \frac{2}{4} \sqrt{5} \right) - \left( \frac{6}{4} - \frac{2}{4} \sqrt{5} \right) \right\} & = & \frac{1}{\sqrt{5}} ・ \frac{4}{4} \sqrt{5} & = & 1 \\ & F_{3} = & \frac{1}{\sqrt{5}} \left\{ \left( \frac{16}{8} + \frac{8}{8} \sqrt{5} \right) - \left( \frac{16}{8} - \frac{8}{8} \sqrt{5} \right) \right\} & = & \frac{1}{\sqrt{5}} ・ \frac{16}{8} \sqrt{5} & = & 2 \\ & F_{4} = & \frac{1}{\sqrt{5}} \left\{ \left( \frac{56}{16} + \frac{24}{16} \sqrt{5} \right) - \left( \frac{56}{16} - \frac{24}{16} \sqrt{5} \right) \right\} & = & \frac{1}{\sqrt{5}} ・ \frac{48}{16} \sqrt{5} & = & 3 \\ & F_{5} = & \frac{1}{\sqrt{5}} \left\{ \left( \frac{176}{32} + \frac{80}{32} \sqrt{5} \right) - \left( \frac{176}{32} - \frac{80}{32} \sqrt{5} \right) \right\} & = & \frac{1}{\sqrt{5}} ・ \frac{160}{32} \sqrt{5} & = & 5 \end{eqnarray*}

こうして見ていくと、二項定理からも明らかなように \(( 1 / 2 + \sqrt{5} / 2 )^{n} - ( 1 / 2 - \sqrt{5} / 2 )^{n}\) を計算すると、有理数は打ち消しあって無理数( \(\sqrt{5}\) を含む項 )だけが残ります。
従って、\(( \frac{1}{2} + \frac{\sqrt{5}}{2} ) ^ {n}\) の \(\sqrt{5}\) の係数を2倍すると \(F_{n}\) と等しくなります。

では、どうしたら \(( \frac{1}{2} + \frac{\sqrt{5}}{2} ) ^ {n}\) の \(\sqrt{5}\) の係数を求められるのでしょうか?

\(「(有理数 + 有理数・\sqrt{5})^{n}」\) は必ず \(「有理数 + 有理数・\sqrt{5}」\) の形になるので、 こちらの記事 にも書いてあるように、虚数の計算のように「有理数」と \(「有理数・\sqrt{5}」\) をわけて計算すれば良さそうです。 例えば \(「5 + 3\sqrt{5}」\) を (5, 3) と表現することにします。


《 下準備 》

ここで、ちょっと寄り道をします。

SICP の中で、 bn を高速に計算する方法として、次のような手続きが紹介されています。

(define (fast-expt b n)
  (cond ((= n 0) 1)
        ((even? n) (square (fast-expt b (/ n 2))))
        (else (* b (fast-expt b (- n 1))))))

でも、この方法は「乗算」だけに限定する必要はないんです。結合則を満す演算なら何でも使えそうです。
そこで『素因数分解と素数判定』という本に載っていたアルゴリズムを参考に次のような関数を作ってみました。

-----------------------------------------
-- 冪乗法
-----------------------------------------
-- e の初期値 : 関数 f の単位元
-- ex : 2^100 == power (*) 2 100 1
--
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)

この関数を使うと \(「a^{n}」\) だけでなく、 素数判定などにでてくる \(「a^{n} \pmod{m}」\) のような計算も高速に行えます。

--
-- a^n (mod m) を計算する
--
powerMod :: (Integral a, Integral t) => t -> a -> t -> t
powerMod a n m = power (\x y -> mod (x * y) m) a n 1




《 一般項からフィボナッチ数を求める関数 》

話をフィボナッチに戻します。

\(a_{1}, a_{2}, b_{1}, b_{2}\) を有理数とすると、

\[( a_{1} + b_{1} \sqrt{5} ) ( a_{2} + b_{2} \sqrt{5} ) = ( a_{1}a_{2} + 5b_{1}b_{2} ) + ( a_{1}b_{2} + a_{2}b_{1} ) \sqrt{5} \]

となることから、次のような関数を定義します。

f (a1, b1) (a2, b2) = (a1 * a2 + 5 * b1 * b2, a1 * b2 + a2 * b1)

この関数と前述の「power 関数」を使えば、次のような「一般項からフィボナッチ数を求める関数」を定義することができます(有理数を扱うので「import Data.Ratio」が必要です)。

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)

fibonacci :: Int -> Integer
fibonacci 0 = 0
fibonacci 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)

「power 関数」がかなり強力なので、この「fibonacci 関数」も高速です。

GHCi 上で実行したら、次のような結果になりました。 (実際に計算させると 20898764 桁の数になるので、結果の表示に時間がかからないように "mod" を付け加えました)。

> mod (fib9 100000000) 3
0
(12.32 secs, 1175839488 bytes)

100,000,000 番目のフィボナッチ数が 12 秒ほどで計算できるのだから、なかなかのものではないでしょうか。

« 2012年10月 | トップページ | 2016年7月 »

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            
フォト

最近のトラックバック

無料ブログはココログ