フィボナッチ数列の一般項を計算する(※ただし有理数に限る)

フィボナッチ数列という有名な数列がある。 1 1 2 3 5 8 … と続くアレだ。漸化式でいうと、

\begin{eqnarray*} F_1 &=& 1 \\ F_2 &=& 1 \\ F_{n+2} &=& F_n + F_{n+1} \ \ (1 \leq n) \end{eqnarray*}

っていうこいつである。

で、こいつの一般項というのが知られている。それは以下の通りである。…らしい。

F_n = \(frac{1}{\sqrt{5}}\biggl\{\biggl(\frac{1 + \sqrt{5}}{2}\biggr)^n – \biggl(\frac{1 – \sqrt{5}}{2}\biggr)^n\biggr\}\)

え、えー？ いやちょっと待ってよだってほら整数の数列だし一般項に無理数入ってるとか大分意味わかんないんですけどこれ大丈夫なの本当に？

うーん。確認してみよう。

fibDouble :: Integer -> Double
fibDouble n = ( ((1 + sqrt 5)/2)^n - ((1 - sqrt 5)/2)^n ) / sqrt 5

main :: IO ()
main = mapM_ (print . fibDouble) [1..10]

runghc すると：

1.0
1.0
2.0
3.0
5.0
8.0
13.0
21.0
34.0
54.99999999999999

うん。あって…る。微妙に誤差があるのだがまあ四捨五入すれば十分使えるレベルのハズ。誤差が怪しかったひとつ前の 9 番目からを改めてやってみよう。

fibDouble :: Integer -> Integer
fibDouble n = round $ ( ((1 + sqrt 5)/2)^n - ((1 - sqrt 5)/2)^n ) / sqrt 5

main :: IO ()
main = mapM_ (print . fibDouble) [9..15]

えいや

34
55
89
144
233
377
610

良いのでは？

一般項を使っているので、 n 番目一個だけを知りたい、という要求においてはこの実装は確かにそれなりに速い。ただ内部の計算が浮動小数なので、100 番目くらいを求めようとするともうもはや求まらない。

main = print $ fibDouble 100

実行すると

354224848179261800448

あってる？　いいえ、残念。100 番目の本当のフィボナッチ数は 354224848179261915075 なのだ。(※ちなみにこの数を Google さんで検索すると、本記事執筆時点で 7,000 件以上ひっかかってそれなりに有名な数であることがわかる。)

もう浮動小数なんかダメだ。信用ならない。まったく、たかが 100 番目くらいで正しい答えを得られないようでは”なってない”といわざるを得ないではないか。やっぱりここは誤差なしで扱えるイケメン、じゃなかった、有理数に限る。

でも、どうやって？ 元の式に無理数が含まれている以上、有理数だけで計算するなんてできっこないじゃないか。

計算方針を考える

落ち着こう。確かに無理数は避けて通れない。しかし、我々は無理数どころか、”もっとスゴい”数を扱ったことがあるはずだ。そう、二乗すると -1 になっちゃう愛に溢れた”アレ”を含む数、複素数である。

複素数をどう扱ったか思い出そう。次のように表現したはずだ。

$$ a + bi $$

虚数単位 \( i \) はもはや無理数ですらない。それでも我々は扱ってきた。ならばたかが \( \sqrt{5} \) 程度が扱えないわけがない。

というわけで、複素数の \( i \) を \( \sqrt{5} \) に置き換えた形で今回の式を取り扱うことにしよう。

データ構造を作る

なんとなく Haskell を使っているので真面目にやるのなら専用の型を用意する場面なのかもしれないが、まぁ、いいや、めんどくせぇ。タプルで十分だ。

import Data.Ratio

type FibNum = (Rational, Rational)

そうそう、Haskell では標準で Rational という有理数を扱う型が用意されている。分母分子を区切るために % 記号を用いることで簡単に有理数を扱える大変イケてるヤツだ。1 % 4 と書けば \( \frac{1}{4} \) を意味する。Rational の分母分子は Integer なので、桁数についても計算資源が許す限りはいくらでも扱える。

さて、この FibNum こと Rational の二要素からなるタプルは、左に \( \sqrt{5} \) が付かない項を、右に \( \sqrt{5} \) が付く項を格納することしよう。つまり (1, 1) と書けば \( 1 + \sqrt{5} \) のこと。(0, 1) と書けば \( \sqrt{5} \) のこと。(1 % 2, 1 % 3) と書けば \( \frac{1}{2} + \frac{1}{3} \sqrt{5} \) のことを表す。

式を書き下す

そんじゃ、素直に一般項を表現してみよう。一般項はこうだった。

$$ F_n = \frac{1}{\sqrt{5}}\biggl\{\biggl(\frac{1 + \sqrt{5}}{2}\biggr)^n – \biggl(\frac{1 – \sqrt{5}}{2}\biggr)^n\biggr\} $$

Haskell 語ではこんな感じか。

fibMain :: Integer -> FibNum
fibMain n = (
    (((1,  1) `fibDiv` (2, 0)) `fibExp` n)
    `fibMin`
    (((1, -1) `fibDiv` (2, 0)) `fibExp` n)
    ) `fibDiv` (0, 1)

OK。さすがにちょっと見づらいがまあ仕方ない。でも fibDiv として素直に除算を除算のまま書き下してしまった。除算は \( 0 \) で割れないとか面倒なこともあるので、乗算の形にしておきたい。

まず、(1, 1) `fibDiv` (2, 0) は要するに \( \frac{1 + \sqrt{5}}{2} \) のことだが、こんなものは \( \frac{1}{2} + \frac{1}{2}\sqrt{5} \)、つまり (1 % 2, 1 % 2) としてしまえば良い。

後ろ側の \( \sqrt{5} \) で割る処理は、逆数であるところの \( \frac{1}{\sqrt{5}} \) 、つまり \( \frac{\sqrt{5}}{5} \) を掛ければ良い。\( \frac{\sqrt{5}}{5} \) ってことは \( 0 + \frac{1}{5}\sqrt{5} \) だから、ここでの表現では (0, 1 % 5) ってことだ。

書き直そう。

fibMain :: Integer -> FibNum
fibMain n = (fibExp (1 % 2, 1 % 2) n `fibMin` fibExp (1 % 2, -1 % 2) n) `fibMul` (0, 1 % 5)

これで除算がなくなった。減算は難しくないのでそのまま実装してしまおう。

演算子を実装する : 累乗

まずやらなきゃならないのは累乗だが、これは乗算の繰り返しの形で表現できる。素直に書くとこうだが：

fibExp :: FibNum -> Integer -> FibNum
fibExp _ 0 = (1, 0)
fibExp a 1 = a
fibExp a n = fibMul a (fibExp a (n - 1))

\( n \) が大きいとそのまま \(n – 1\) 回の乗算をすることになってちょっとばかり遅い。二乗の結果が使えるところはどんどんそれを使って計算させることにしよう。乗算の回数が最大でも \( 2 \log{n} \) 回で済む。

fibExp a n
    | even n    = ex
    | otherwise = ex `fibMul` a
    where ex = fibExp (fibMul a a) (n `div` 2)

演算子を実装する : 乗算

FibNum の乗算とは何かと言うと、\( (a + b\sqrt{5})(c + d\sqrt{5}) \) ってことで、つまり、

\begin{eqnarray*} & & (a + b\sqrt{5})(c + d\sqrt{5}) \\ &=& ac + ad\sqrt{5} + bc\sqrt{5} + 5bd \\ &=& (ac + 5bd) + (ad + bc)\sqrt{5} \end{eqnarray*}

Haskell 語で書くと：

fibMul :: FibNum -> FibNum -> FibNum
fibMul (a, b) (c, d) = (a*c + 5*b*d, a*d + b*c)

演算子を実装する : 減算

今回何故か加算ではなく減算だが、まぁいいだろう。目的の式がたまたまそうだっただけだ。これは別に難しくない：

fibMin :: FibNum -> FibNum -> FibNum
fibMin (a, b) (c, d) = (a - c, b - d)

同類項それぞれで計算するだけだ。

使ってみよう！

ここまでくればあとは最初に定義した fibMain を呼ぶだけだ。

main :: IO ()
main = mapM_ (print . fibMain) [1..10]

実行すると：

(1 % 1,0 % 1)
(1 % 1,0 % 1)
(2 % 1,0 % 1)
(3 % 1,0 % 1)
(5 % 1,0 % 1)
(8 % 1,0 % 1)
(13 % 1,0 % 1)
(21 % 1,0 % 1)
(34 % 1,0 % 1)
(55 % 1,0 % 1)

おおっと…。表記がちょっと見づらい。

フィボナッチ数列は元々整数の数列なので、当然計算結果に無理数であるところの \( \sqrt{5} \) は絶対に含まれない。だから、 \( a + b \sqrt{5} \) の \( b \) は絶対に \( 0 \) になるのだ。それに、整数ということは有理数 \( a \) の分母も絶対に \( 1 \) になる。余分なコトを出力すると見づらいので、その辺は削ってしまおう。

fib :: Integer -> Integer
fib n = numerator a
    where (a, _) = fibMain n

main も呼ぶ関数を直して

main :: IO ()
main = mapM_ (print . fib) [1..10]

もっかい

1
1
2
3
5
8
13
21
34
55

どや！

イケメンの実力を確認する

ここまでくれば 100 番目だって正確に求められる。そう、有理数ならね。

main = print $ fib 100

ほれ。

354224848179261915075

別に 10,000 番目だっていける。

33644764876431783266621612005107543310302148460680063906564769974680081442166662368155595513633734025582065332680836159373734790483865268263040892463056431887354544369559827491606602099884183933864652731300088830269235673613135117579297437854413752130520504347701602264758318906527890855154366159582987279682987510631200575428783453215515103870818298969791613127856265033195487140214287532698187962046936097879900350962302291026368131493195275630227837628441540360584402572114334961180023091208287046088923962328835461505776583271252546093591128203925285393434620904245248929403901706233888991085841065183173360437470737908552631764325733993712871937587746897479926305837065742830161637408969178426378624212835258112820516370298089332099905707920064367426202389783111470054074998459250360633560933883831923386783056136435351892133279732908133732642652633989763922723407882928177953580570993691049175470808931841056146322338217465637321248226383092103297701648054726243842374862411453093812206564914032751086643394517512161526545361333111314042436854805106765843493523836959653428071768775328348234345557366719731392746273629108210679280784718035329131176778924659089938635459327894523777674406192240337638674004021330343297496902028328145933418826817683893072003634795623117103101291953169794607632737589253530772552375943788434504067715555779056450443016640119462580972216729758615026968443146952034614932291105970676243268515992834709891284706740862008587135016260312071903172086094081298321581077282076353186624611278245537208532365305775956430072517744315051539600905168603220349163222640885248852433158051534849622434848299380905070483482449327453732624567755879089187190803662058009594743150052402532709746995318770724376825907419939632265984147498193609285223945039707165443156421328157688908058783183404917434556270520223564846495196112460268313970975069382648706613264507665074611512677522748621598642530711298441182622661057163515069260029861704945425047491378115154139941550671256271197133252763631939606902895650288268608362241082050562430701794976171121233066073310059947366875

って、でかっ。(※右がはみ出て切れているかもしれないが、2,090 桁の数である)

でもっ、さらにっ！ 1,000,000 番目だっていける！

1953282128 ...(中略)... 8242546875

…のだが、実に 208,989 桁に及ぶ数が出力されるので全桁載せるのはやめておこう。どうせ表示されない。

完成品

ソースコード全体はこんな感じ。main はコマンドライン引数にて n 番目を指示させるようにしてある。まぁ数字じゃないコマンドライン引数渡すとそのままエラー吐いて終了する不親切設計だけど、そこは大目に見てもらうってことで。

import System.Environment (getArgs)
import Data.Ratio

type FibNum = (Rational, Rational)

fibMain :: Integer -> FibNum
fibMain n = (fibExp (1 % 2, 1 % 2) n `fibMin` fibExp (1 % 2, -1 % 2) n) `fibMul` (0, 1 % 5)

fibExp :: FibNum -> Integer -> FibNum
fibExp _ 0 = (1, 0)
fibExp a 1 = a
fibExp a n
    | even n    = ex
    | otherwise = ex `fibMul` a
    where ex = fibExp (fibMul a a) (n `div` 2)

fibMul :: FibNum -> FibNum -> FibNum
fibMul (a, b) (c, d) = (a*c + 5*b*d, a*d + b*c)

fibMin :: FibNum -> FibNum -> FibNum
fibMin (a, b) (c, d) = (a - c, b - d)

fib :: Integer -> Integer
fib n = numerator a
    where (a, _) = fibMain n

main :: IO ()
main = mapM_ (print . fib . read) =<< getArgs

で、実は

この実装、さすがに一般項を使っているだけあって n 番目を一発で求める分にはまあそこそこ速いのだけど、それでももっと速い実装が既に存在する。有名な数列なだけあって、行列を用いて漸化式を変形させた方法等があるようだ。ま、正確かつそれなりに高速に計算できたってことで、個人的にはそれなりに満足である。ちなみに手元の環境では 10,000,000 番目を求めるのに 4 秒程度、100,000,000 番目で 1.6 分程度であった。