LoginSignup
156
95

More than 5 years have passed since last update.

@mod_poppo

Haskellでフィボナッチ数列 〜Haskellで非実用的なコードを書いて悦に入るのはやめろ〜

Last updated at Posted at 2019-01-19

フィボナッチ数列の定義と素朴な実装

フィボナッチ数列は漸化式$$
F_0=0,\quad
F_1=1,\quad
F_n=F_{n-2}+F_{n-1}\ (n\ge 2)
$$で定まる数列で、最初の方をいくつか挙げると 0, 1, 1, 2, 3, 5, 8, … となります。（流儀によっては1から始めることもあります）

では、Haskellでフィボナッチ数を計算するコードを書いてみましょう。

VerySlowFib.hs
-- 素朴なコード
fib :: Int -> Integer
fib 0 = 0
fib 1 = 1
fib n = fib (n - 2) + fib (n - 1)

数学的な定義をそのまま書き下した、Haskellらしい美しいコードですね！:innocent:

……なーんて感想を抱いた人はHaskell初心者です。計算機で動かすコードを「数学的な美しさ」で語ってどうするんですか？？？

美しいバラにはトゲがあるように、美しいHaskellコードには罠があります。

この「素朴なコード」は非常に遅く、実用に耐えません。（この素朴なコードの出番があるとしたら「遅いコードの例」としてでしょう）

素朴なコードの問題点

例として、「素朴なコード」で fib 4 を計算してみましょう。

fib 4 = fib 2 + fib 3
      = (fib 0 + fib 1) + fib 3
      = (0 + fib 1) + fib 3
      = (0 + 1) + fib 3
      = 1 + fib 3
      = 1 + (fib 1 + fib 2)
      = 1 + (1 + fib 2)
      = 1 + (1 + (fib 0 + fib 1))
      = 1 + (1 + (0 + fib 1))
      = 1 + (1 + 1)
      = 1 + 2
      = 3

やたら手数を食っている（fib 2 を2回計算している）のがわかりますか？一度計算した値を覚えておけばもっとシュッと計算できそうです。つまり、

fib 2 = fib 0 + fib 1
      = 0 + fib 1
      = 0 + 1
      = 1
fib 3 = fib 1 + fib 2
      = 1 + 1
      = 2
fib 4 = fib 2 + fib 3
      = 1 + 2
      = 3

とすれば fib 2 を1回しか計算しなくて済みますね。「1回計算した値を覚えておく」のにはメモ化と呼ばれるテクニックが使われます。

Haskellコンパイラー (GHC) は魔法じゃないので関数が純粋だからといって自動的にメモ化したりはしてくれません！

リストを使ったメモ化

その辺のプログラミング言語でメモ化と言ったらハッシュテーブルとかを使うと思いますが、Haskellではハッシュテーブルのような可変なデータ構造を手軽に扱うことはできません（ライブラリーを使ったりSTモナドを使ったりする必要がある）。ここでは、遅延リストを使った簡易的なメモ化を実装してみます。

ListMemo.hs
-- リストを使ったメモ化
fib :: Int -> Integer
fib 0 = 0
fib 1 = 1
fib n = fibList !! (n - 2) + fibList !! (n - 1)

fibList :: [Integer]
fibList = map fib [0..]

fib n の右辺で fib 自身を直接呼び出すのではなく、 fibList を参照していることに注意してください。fibList の定義から fibList !! n == fib n ですが、 fibList の方は一度計算した値を覚えています。これならより高速にフィボナッチ数を計算できそうです。

この方法の欠点は、リストに対して !! でアクセスしていることです（連結リストへの添字でのアクセスは添字に比例した時間がかかる）。メモ化の方法として二分木やハッシュテーブルを使えばもう少しマシになるかもしれませんが、そもそも、計算の際に直前2項しか使わないのに汎用的なメモ化を行うのは無駄ですね。

zipWith を使った定義（再帰的なリスト定義）

Haskellでフィボナッチ数と言えば、 zipWith を使った

ListFib.hs
module ListFib where

-- zipWithを使った版
fib :: [Integer]
fib = 0 : 1 : zipWith (+) fib (tail fib)

が有名かと思います。詳しい解説はしません1が、この無限リストを評価すると [0,1,1,2,3,5,8,...] という風にフィボナッチ数になります。

しかし、このコードには遅延評価絡みの問題があり、次のように改変すると高速化します：

ListFib.hs
-- リストのhead側を正格評価する版
zipWith' :: (a -> b -> c) -> [a] -> [b] -> [c]
zipWith' f (x:xs) (y:ys) = let z = f x y
                           in z `seq` (z : zipWith' f xs ys)

-- zipWithを使った版（正格評価）
fib' :: [Integer]
fib' = 0 : 1 : zipWith' (+) fib' (tail fib')

詳細については、すでにいくつか書かれている記事を参照してください：

再帰関数で計算する

さっきの「zipWithを使った定義」とほぼ同等の処理を（リストを使わず）再帰関数で書くと次のようになります：

IterFib.hs
{-# LANGUAGE BangPatterns #-}
module IterFib where

-- 再帰関数で計算する（遅延評価）
fibL :: Int -> Integer
fibL i = loop i 0 1
  where
    loop 0 a b = a
    loop i a b = loop (i - 1) b (a + b)

-- 再帰関数で計算する
fib :: Int -> Integer
fib i = loop i 0 1
  where
    loop 0 a !b = a
    loop i a b = loop (i - 1) b (a + b)

2通り書きましたが、GHCの最適化を有効にした場合は IterFib.fibL よりも IterFib.fib の方が高速となります。（IterFib.fibL では loop の b が正格にならず、サンクが作られる）

この「再帰関数で計算する」版 IterFib.fib でn番目のフィボナッチ数を計算するのと「zipWithを使った版（正格評価）」 IterList.fib' でn番目を計算するのとでは、実行時間がほぼ同じとなります。

n番目のフィボナッチ数だけをピンポイントで計算する

「n番目までのフィボナッチ数を全て計算する」場合はさっきの「zipWithを使った版（正格評価）」と「再帰関数で計算する版」で話は終わりですが、「n番目のフィボナッチ数だけをピンポイントで計算する（n-1番目までのフィボナッチ数は必要ない）」場合はさらに早い方法があります。

（逆に、これ以降で述べる方法はn番目までのフィボナッチ数を全て計算する場合には適さないので注意してください。適材適所です。）

フィボナッチ数列の一般項（ビネの公式）を使う

フィボナッチ数列には一般項があります（参考：Wikipediaの項目）。この公式はビネ (Binet) の公式と呼ばれているようです。
$$
F_n=\frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^n-\left(\frac{1-\sqrt{5}}{2}\right)^n\right)
$$

これを使うとn番目のフィボナッチ数を一発で計算できます！　……机上の、数学的な議論ならここで終わっても良さそうですが、あいにく我々は計算機上で実装しないと気が済まない種族ですので、話はまだ終わりません。

ビネの公式では $\sqrt{5}$ が使われていますが、計算機でルート5という数はどうやって表現したものでしょうか。ここでは何通りか検討してみましょう。

浮動小数点数を使う

まずは倍精度浮動小数点数を使ってビネの公式を計算してみましょう。

BinetDouble.hs
module BinetDouble where

phi, psi :: Double
phi = (1 + sqrt 5) / 2
psi = (1 - sqrt 5) / 2

-- ビネの公式を浮動小数点数で・その1
fib :: Int -> Integer
fib i = round $ (phi^i - psi^i) / sqrt 5

-- ビネの公式を浮動小数点数で・その2
fibP :: Int -> Integer
fibP i = round $ (phi ** fromIntegral i - psi ** fromIntegral i) / sqrt 5

その1（BinetDouble.fib）は浮動小数点数のべき乗にHaskellの (^) 演算子を使っています。(^) は浮動小数点数の乗算を繰り返し行うことでべき乗を計算します。一方、その2（BinetDouble.fibP）は (**) 演算子を使ってべき乗を計算します。こっちは他の言語では pow 関数みたいな名前がついているやつです。

さて、浮動小数点数を使う欠点ですが、有効数字が限られているので、nが大きくなると正しい結果が得られなくなります。

具体的には、倍精度の場合は BinetDouble.fibP は n=71 で、 BinetDouble.fib は n=76 で破綻（正しい整数値が得られない）します。この辺のフィボナッチ数は $F_{71} = 308061521170129$ というような10進で15桁の数になるので、Doubleの精度が15桁か16桁であることを考えれば自然ですね。

というわけで浮動小数点数を使ってフィボナッチ数列の一般項を計算する方法は使い物になりません。

なになに、「俺は n≤70 しか使わないから困らないぞ」ですって？　言いづらいのですが、そもそも n≤100 しか必要ないのだったら、実行時に計算しなくても、値を全部ソースコードに埋め込んでおけばいいのではないのでしょうか。

代数拡大体 Q(√5) を使う

ルートを含む数を計算機で表す方法は浮動小数点数だけではありません。筆者の書いた「週刊 代数的実数を作る」ではルートを含む数、もっと一般に代数的数と呼ばれるクラスの数を計算機上で正確に取り扱う方法を解説しています。気になった方は読んでみてください。サンプルコードはHaskellで書かれています。

閑話休題。

今回は $\sqrt{5}$ さえ扱えれば良いので、「週刊 代数的実数を作る」に書いたような一般の代数的数を使うのはオーバーキルです。ここは、 $\mathbf{Q}$ の2次拡大$$
\mathbf{Q}(\sqrt{5})=\{a+b\sqrt{5}\mid a,b\in\mathbf{Q}\}
$$を使っておけば良いでしょう。実装例は次のようになります：

BinetSqrt5.hs
module BinetSqrt5 where
import Data.Ratio

-- Ext_sqrt5 a b = a + b * sqrt 5
data Ext_sqrt5 a = Ext_sqrt5 !a !a deriving (Eq,Show)

instance Num a => Num (Ext_sqrt5 a) where
  Ext_sqrt5 a b + Ext_sqrt5 a' b' = Ext_sqrt5 (a + a') (b + b')
  Ext_sqrt5 a b - Ext_sqrt5 a' b' = Ext_sqrt5 (a - a') (b - b')
  negate (Ext_sqrt5 a b) = Ext_sqrt5 (negate a) (negate b)
  Ext_sqrt5 a b * Ext_sqrt5 a' b' = Ext_sqrt5 (a * a' + 5 * b * b') (a * b' + b * a')
  fromInteger n = Ext_sqrt5 (fromInteger n) 0

instance Fractional a => Fractional (Ext_sqrt5 a) where
  recip (Ext_sqrt5 a b) = let s = a * a - 5 * b * b
                          in Ext_sqrt5 (a / s) (- b / s)
  fromRational x = Ext_sqrt5 (fromRational x) 0

phi, psi, sqrt5 :: Ext_sqrt5 Rational
phi = Ext_sqrt5 (1/2) (1/2)
psi = Ext_sqrt5 (1/2) (-1/2)
sqrt5 = Ext_sqrt5 0 1

fib :: Int -> Integer
fib i = case (phi^i - psi^i) / sqrt5 of
          Ext_sqrt5 x 0 | denominator x == 1 -> numerator x
          x -> error $ "calculation error: fib " ++ show i ++ " = " ++ show x

Ext_sqrt5 Rational が $\mathbf{Q}(\sqrt{5})$ に相当する型です。これを使うと、黄金比 $\phi=\frac{1+\sqrt{5}}{2}$ は Ext_sqrt5 (1/2) (1/2) と表現できます。

このコードでは、 phi^i および psi^i の部分（べき乗計算）で主に時間がかかっています。浮動小数点数の場合はべき乗計算は一瞬だったかもしれませんが、多倍長整数・有理数の組で正確な計算をするとべき乗の部分にも一定のコストがかかります。

代数拡大体 Q((1+√5)/2) を使う

さっきは有理数体に √5 を添加しましたが、代わりに $\phi=\frac{1+\sqrt{5}}{2}$ を添加する方法も考えてみます。$\phi$ を使うとビネの公式は
$$
F_n=\frac{\phi^n-(1-\phi)^n}{2\phi-1}
$$
となり、分子の部分が整数演算のみで済みそうに思えます。実装は次のようになります：

BinetPhi.hs
module BinetPhi where
import Data.Ratio

{-
Ext_phi a b = a + b * phi
where phi is a root of
  phi^2 - phi - 1 = 0
-}
data Ext_phi a = Ext_phi !a !a deriving (Eq,Show)

instance Num a => Num (Ext_phi a) where
  Ext_phi a b + Ext_phi a' b' = Ext_phi (a + a') (b + b')
  Ext_phi a b - Ext_phi a' b' = Ext_phi (a - a') (b - b')
  negate (Ext_phi a b) = Ext_phi (negate a) (negate b)
  Ext_phi a b * Ext_phi a' b' = let bb' = b * b'
                                in Ext_phi (a * a' + bb') (a * b' + b * a' + bb')
  fromInteger n = Ext_phi (fromInteger n) 0

instance Fractional a => Fractional (Ext_phi a) where
  recip (Ext_phi a b) = let s = a * a + a * b - b * b
                        in Ext_phi ((a + b) / s) (- b / s)
  fromRational x = Ext_phi (fromRational x) 0

phi, psi :: (Num a) => Ext_phi a
phi = Ext_phi 0 1
psi = 1 - phi

fib :: Int -> Integer
fib i = case (phi^i - psi^i) / (phi - psi) of
          Ext_phi x 0 | denominator x == 1 -> numerator x
          x -> error $ "calculation error: fib " ++ show i ++ " = " ++ show x

fibI :: Int -> Integer
fibI i = case phi^i - psi^i of
          Ext_phi mx y | 2 * (- mx) == y -> - mx
          x -> error $ "calculation error: fib " ++ show i ++ " * sqrt 5 = " ++ show x

BinetPhi.fib は phi^i - psi^i の部分を有理数係数 $\mathbf{Q}(\frac{1+\sqrt{5}}{2})$ で、 BinetPhi.fibI は整数係数 $\mathbf{Z}[\frac{1+\sqrt{5}}{2}]$ で計算しています。

実行時間ですが、 BinetPhi.fib および BinetPhi.fibI ともに、 BinetSqrt5.fib で計算したのとほとんど同じでした（後述）。整数係数になるとはいえ、 $\frac{1+\sqrt{5}}{2}$ を添加する方は乗算が複雑になっていることが関係しているのでしょうか。

行列のべき乗を使う

フィボナッチ数の漸化式は、行列によって次のように書けます：
$$
\begin{pmatrix}
F_{n}\\F_{n+1}
\end{pmatrix}=
\begin{pmatrix}
0&1\\
1&1
\end{pmatrix}
\begin{pmatrix}
F_{n-1}\\F_{n}
\end{pmatrix}
$$

この関係式を繰り返し使うと、n番目およびn+1番目のフィボナッチ数は行列のべき乗を使って次のように計算できます：
$$
\begin{pmatrix}
F_{n}\\F_{n+1}
\end{pmatrix}=
\begin{pmatrix}
0&1\\
1&1
\end{pmatrix}^n
\begin{pmatrix}
F_{0}\\F_{1}
\end{pmatrix}=
\begin{pmatrix}
0&1\\
1&1
\end{pmatrix}^n
\begin{pmatrix}
0\\1
\end{pmatrix}
$$

つまり、 $F_n$ は行列 $\begin{pmatrix}
0&1\\
1&1
\end{pmatrix}^n$ の右上の成分、ということですね。行列のべき乗は、 $n$ の部分を二進展開する計算法で、そこそこの速度で計算できます2。

実装例は次のようになります：

{-# LANGUAGE BangPatterns #-}
module MatFib where

{-
Mat2x2 a b c d
  = / a b \
    \ c d /
-}
data Mat2x2 a = Mat2x2 !a !a !a !a

matMul :: (Num a) => Mat2x2 a -> Mat2x2 a -> Mat2x2 a
matMul (Mat2x2 a b c d) (Mat2x2 a' b' c' d')
  = Mat2x2 (a * a' + b * c') (a * b' + b * d')
           (c * a' + d * c') (c * b' + d * d')

matPow :: (Num a) => Mat2x2 a -> Int -> Mat2x2 a
matPow _ 0 = Mat2x2 1 0 0 1
matPow m i = loop m m (i - 1)
  where
    loop acc !_ 0 = acc
    loop acc m 1 = acc `matMul` m
    loop acc m i = case i `quotRem` 2 of
               (j,0) -> loop acc (m `matMul` m) j
               (j,_) -> loop (acc `matMul` m) (m `matMul` m) j

fib :: Int -> Integer
fib i = let Mat2x2 a b c d = matPow (Mat2x2 0 1 1 1) i
        in b

この MatFib.fib と先ほどの BinetSqrt5.fib, BinetPhi.fib, BinetPhi.fibI はほぼ同じ実行時間となりました。

実験

これまでに説明したアルゴリズムを実際に動かして、実行時間を比べてみましょう。

コードは

に置いておきます。ビルド手順は以下です：

$ git clone https://github.com/minoki/fibonacci-hs.git
$ cd fibonacci-hs
$ stack build

stack exec で実行するコマンドの実行時間を測る際のコツですが、 time stack exec ではなく stack exec time とします（筆者の環境では、stack exec自体の実行時間が0.2秒ほどかかるようなので）。

以下、筆者の環境 (MacBook Pro Late 2013, GHC 8.6.3) で実行した結果（実行時間）を貼っておきます。

順番に計算する場合

まずは、最初の40個のフィボナッチ数を列挙してみましょう。「素朴な方法」で50個列挙しようとしたところ、せっかちな筆者には待ちきれないほど時間がかかったので40個にとどめています。

$ # 素朴な方法で n=40 まで計算する。4秒ほど
$ stack exec time fibonacci-hs-exe VerySlow 40
F[0] = 0
F[1] = 1
F[2] = 1
（略）
F[38] = 39088169
F[39] = 63245986
        4.09 real         3.90 user         0.04 sys
$ # リストによるメモ化を適用して n=40 まで計算する。一瞬
$ stack exec time fibonacci-hs-exe ListMemo 40
F[0] = 0
F[1] = 1
F[2] = 1
（略）
F[38] = 39088169
F[39] = 63245986
        0.02 real         0.00 user         0.00 sys

素朴な方法 (VerySlowFib) と リストによるメモ化 (ListMemo) では後者の圧勝です。素朴な方法が実用的じゃないことが改めてわかりました。

次は、リストによるメモ化 (ListMemo) と他のいくつかの方法を比べてみましょう。「フィボナッチ数を10000で割った余り」の最初の20000（にまん）項の和を10000で割った余り、を計算させます（要はフィボナッチ数を最初から順に20000項計算しています）。

$ # リストによるメモ化の場合
$ stack exec time fibonacci-hs-exe ListMemo sum 20000
625
       23.55 real        18.31 user         0.24 sys
$ # zipWith（遅延評価）で定義した場合
$ stack exec time fibonacci-hs-exe List sum 20000
625
        0.03 real         0.01 user         0.00 sys
$ # 毎回再帰で計算する
$ stack exec time fibonacci-hs-exe Iter sum 20000
625
       15.63 real        15.16 user         0.30 sys
$ # Q(√5)でビネの公式を使う
$ stack exec time fibonacci-hs-exe BinetSqrt5 sum 20000
625
        2.00 real         1.97 user         0.01 sys
$ # 行列のべき乗を使う
$ stack exec time fibonacci-hs-exe Mat sum 20000
625
        1.12 real         1.07 user         0.01 sys

zipWith（遅延評価）(ListFib.fib) が一瞬なのに対し、リストによるメモ化 (ListMemo) は20秒以上かかっています。また、ビネの公式や行列のべき乗を使う方法は、1秒から2秒程度と、zipWith（遅延評価）(ListFib.fib) よりも遅くなっています。

今度は、zipWith（遅延評価）(ListFib.fib) とその正格評価版を比べてみましょう。 最初の2000000（にひゃくまん）項を計算させます。

$ # zipWith（遅延評価） (ListFib.fib)
$ stack exec time fibonacci-hs-exe List sum 2000000
625
       49.69 real        46.28 user         0.81 sys
$ # zipWith'（正格評価） (ListFib.fib')
$ stack exec time fibonacci-hs-exe ListS sum 2000000
625
       43.81 real        42.44 user         0.64 sys

若干、正格評価する方が早くなっていますが、そこまで大きな違いではありません。後述の「1000000（ひゃくまん）番目のフィボナッチ数をピンポイントで計算する」実験で5倍の差が出ているのとは対照的です。

ピンポイントに計算する場合

今度は、大きなフィボナッチ数をピンポイントで計算してみます。

50000（ごまん）番目のフィボナッチ数をそれぞれの方法で計算してみましょう。計算結果はものすごく大きくなるので、 $F_{50000}$ そのものではなく10000で割ったあまりを表示させます。

$ stack exec time fibonacci-hs-exe ListMemo atMod 50000
F[50000] === 3125 mod 10000
       56.78 real        55.11 user         0.51 sys
$ stack exec time fibonacci-hs-exe List atMod 50000
F[50000] === 3125 mod 10000
        0.13 real         0.07 user         0.04 sys
$ stack exec time fibonacci-hs-exe ListS atMod 50000
F[50000] === 3125 mod 10000
        0.05 real         0.02 user         0.01 sys
$ stack exec time fibonacci-hs-exe Iter atMod 50000
F[50000] === 3125 mod 10000
        0.03 real         0.01 user         0.00 sys
$ stack exec time fibonacci-hs-exe IterL atMod 50000
F[50000] === 3125 mod 10000
        0.13 real         0.07 user         0.04 sys
$ stack exec time fibonacci-hs-exe BinetDouble atMod 50000
F[50000] === 7216 mod 10000
        0.02 real         0.00 user         0.00 sys
$ stack exec time fibonacci-hs-exe BinetSqrt5 atMod 50000
F[50000] === 3125 mod 10000
        0.02 real         0.00 user         0.00 sys

ListMemoが圧倒的に時間を食っていることがわかります。また、浮動小数点数を使うBinetDoubleは間違った結果を返しています。それ以外はどれも一瞬なので、ListMemo以外にはもう少し大きいフィボナッチ数を計算させて比べることにします。

1000000（ひゃくまん）番目のフィボナッチ数をそれぞれの方法で計算してみましょう。例によって10000で割ったあまりを表示させます。

$ stack exec time fibonacci-hs-exe List atMod 1000000
F[1000000] === 6875 mod 10000
       32.69 real        32.11 user         0.36 sys
$ stack exec time fibonacci-hs-exe ListS atMod 1000000
F[1000000] === 6875 mod 10000
        6.27 real         6.07 user         0.14 sys
$ stack exec time fibonacci-hs-exe Iter atMod 1000000
F[1000000] === 6875 mod 10000
        6.53 real         6.33 user         0.14 sys
$ stack exec time fibonacci-hs-exe IterL atMod 1000000
F[1000000] === 6875 mod 10000
       33.19 real        32.55 user         0.40 sys
$ stack exec time fibonacci-hs-exe BinetDouble atMod 1000000
F[1000000] === 7216 mod 10000
        0.02 real         0.00 user         0.00 sys
$ stack exec time fibonacci-hs-exe BinetSqrt5 atMod 1000000
F[1000000] === 6875 mod 10000
        0.06 real         0.04 user         0.00 sys
$ stack exec time fibonacci-hs-exe BinetPhi atMod 1000000
F[1000000] === 6875 mod 10000
        0.07 real         0.04 user         0.00 sys
$ stack exec time fibonacci-hs-exe BinetPhiI atMod 1000000
F[1000000] === 6875 mod 10000
        0.07 real         0.04 user         0.00 sys
$ stack exec time fibonacci-hs-exe Mat atMod 1000000
F[1000000] === 6875 mod 10000
        0.06 real         0.04 user         0.00 sys

遅延評価を使う zipWith (List), Iter.fibL (IterL) がやや遅く（30秒程度）、これらの性格評価版である zipWith' (ListS), Iter.fib (Iter) が6秒ちょっとかかっています。一方、べき乗を使う連中（ビネの公式と、行列）はどれも一瞬で計算が終わっています。

最後に、べき乗を使うアルゴリズムを使って100000000（いちおく）番目のフィボナッチ数（の下4桁）を計算してみましょう。

$ stack exec time fibonacci-hs-exe BinetDouble atMod 100000000
F[100000000] === 7216 mod 10000
        0.02 real         0.00 user         0.00 sys
$ stack exec time fibonacci-hs-exe BinetSqrt5 atMod 100000000
F[100000000] === 6875 mod 10000
       10.26 real         9.86 user         0.30 sys
$ stack exec time fibonacci-hs-exe BinetPhi atMod 100000000
F[100000000] === 6875 mod 10000
        9.52 real         9.15 user         0.28 sys
$ stack exec time fibonacci-hs-exe BinetPhiI atMod 100000000
F[100000000] === 6875 mod 10000
        8.93 real         8.54 user         0.27 sys
$ stack exec time fibonacci-hs-exe Mat atMod 100000000
F[100000000] === 6875 mod 10000
        9.75 real         9.37 user         0.29 sys

浮動小数点数を使うBinetDoubleは一瞬で計算が終わっていますが、そもそも答えが間違っているので論外です。他はどれも9秒前後で、有意な差はなさそうです。

まとめ

  • フィボナッチ数を順番に計算する場合には zipWith'（正格評価）による再帰的な定義 あるいは 再帰関数を使うアルゴリズム (IterFib.fib) が実用的です。
{-# LANGUAGE BangPatterns #-}

-- リストのhead側を正格評価する
zipWith' :: (a -> b -> c) -> [a] -> [b] -> [c]
zipWith' f (x:xs) (y:ys) = let z = f x y
                           in z `seq` (z : zipWith' f xs ys)

-- zipWith'を使った版（正格評価）
fibList :: [Integer]
fibList = 0 : 1 : zipWith' (+) fibList (tail fibList)

-- 再帰関数で計算する
fib :: Int -> Integer
fib i = loop i 0 1
  where
    loop 0 a !b = a
    loop i a b = loop (i - 1) b (a + b) -- ここにフィボナッチ数を使った処理を挟む

  • n番目のフィボナッチ数をピンポイントで計算する場合には「代数拡大体上でビネの公式を使う」あるいは「行列のべき乗を使う」方法が実用的です。（フィボナッチ数の実用性についてはこの記事では考えません:wink:）

  • フィボナッチ数を順番に計算する場合とn番目だけをピンポイントで計算する場合では最適なアルゴリズムが違うことに気をつけましょう。

  • Haskellは魔法じゃねえんだから数学的な定義をそのまま書いても実用的なプログラムにはならないぞ:punch:

2019年1月21日 追記：ブログに続編を書きました。この記事に書いたアルゴリズム（行列、ビネの公式）からもっと（定数倍）速くします。

  1. Qiitaにも解説記事があるようです：Haskellのキモいフィボナッチ数列がやっと理解できたからこれでもかという程に細かく説明してみた #Haskell ↩

  2. べきの部分を二進展開する計算法は、演算にかかる時間が対象に依存しない場合は対数時間ですが、今回は行列の中身が多倍長整数という可変な物体なので、対数時間になるとは限りません。 ↩

156
95
8

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
156
95