« セクシー素数 | トップページ | フィボナッチ数・いろいろ »

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 秒ほどで計算できるのだから、なかなかのものではないでしょうか。

« セクシー素数 | トップページ | フィボナッチ数・いろいろ »

Haskell」カテゴリの記事

数学」カテゴリの記事

コメント

コメントを書く

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

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

トラックバック

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

この記事へのトラックバック一覧です: フィボナッチ数を一般項から求めてみる。:

« セクシー素数 | トップページ | フィボナッチ数・いろいろ »

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

最近のトラックバック

無料ブログはココログ