Fibonacci's closed-form expression, the ST monad, and Haskell

Two recent questions about Fibonacci's closed-form expression (here and here) as well as the HaskellWiki's page about the ST monad motivated me to try and compare two ways of calculating Fibonacci numbers.

The first implementation uses the closed-form expression together with rationals as seen in hammar's answer here (where Fib is a datatype abstracting numbers of the form a+b*√5):

fibRational :: Integer -> Integer
fibRational n = divSq5 $ phi^n - (1-phi)^n
  where
    phi = Fib (1/2) (1/2)
    divSq5 (Fib 0 b) = numerator b

The second implementation is from the HaskellWiki's page about the ST monad, with some added strictness that was necessary in order to avoid a stack overflow:

fibST :: Integer -> Integer
fibST n | n < 2 = n
fibST n = runST $ do
    x <- newSTRef 0
    y <- newSTRef 1
    fibST' n x y
  where
    fibST' 0 x _ = readSTRef x
    fibST' !n x y = do
      x' <- readSTRef x
      y' <- readSTRef y
      y' `seq` writeSTRef x y'
      x' `seq` writeSTRef y (x'+y')
      fibST' (n-1) x y

For reference, here's also the full code that I used for testing:

{-# LANGUAGE BangPatterns #-}

import Data.Ratio
import Data.STRef.Strict
import Control.Monad.ST.Strict
import System.Environment

data Fib =
  Fib !Rational !Rational
  deriving (Eq, Show)

instance Num Fib where
  negate (Fib a b) = Fib (-a) (-b)
  (Fib a b) + (Fib c d) = Fib (a+c) (b+d)
  (Fib a b) * (Fib c d) = Fib (a*c+5*b*d) (a*d+b*c)
  fromInteger i = Fib (fromInteger i) 0
  abs = undefined
  signum = undefined

fibRational :: Integer -> Integer
fibRational n = divSq5 $ phi^n - (1-phi)^n
  where
    phi = Fib (1/2) (1/2)
    divSq5 (Fib 0 b) = numerator b

fibST :: Integer -> Integer
fibST n | n < 2 = n
fibST n = runST $ do
    x <- newSTRef 0
    y <- newSTRef 1
    fibST' n x y
  where
    fibST' 0 x _ = readSTRef x
    fibST' !n x y = do
   开发者_开发知识库   x' <- readSTRef x
      y' <- readSTRef y
      y' `seq` writeSTRef x y'
      x' `seq` writeSTRef y (x'+y')
      fibST' (n-1) x y

main = do
  (m:n:_) <- getArgs 
  let n' = read n
      st = fibST n'
      rt = fibRational n'
  case m of
    "st" -> print st
    "rt" -> print rt
    "cm" -> print (st == rt)

Now it turns out that the ST version is significantly slower than the closed-form version, although I'm not a hundred percent sure why:

# time ./fib rt 1000000 >/dev/null
./fib rt 1000000 > /dev/null  0.23s user 0.00s system 99% cpu 0.235 total

# time ./fib st 1000000 >/dev/null
./fib st 1000000 > /dev/null  11.35s user 0.06s system 99% cpu 11.422 total

So my question is: Can someone help me understand why the first implementation is so much faster? Is it algorithmic complexity, overhead or something else entirely? (I checked that both functions yield the same result). Thanks!

You are comparing very different versions here. To make it fair, here is an implementation that is equivalent to the ST solution you give, but in pure Haskell:

fibIt :: Integer -> Integer
fibIt n | n < 2 = n
fibIt n = go 1 1 (n-2)
  where go !_x !y  0 = y
        go !x  !y  i = go y (x+y) (i-1)

This one seems to perform exactly as good or bad as the ST version (both 10s here). The runtime is most likely dominated by all the Integer additions, overhead is therefore too low to be measurable.

First, the two implementations use two very different algorithms with different asymptotic complexity (well, depending on what the complexity of the Integer operations are). Second, the st implementation is using references. References are (comparatively) slow in ghc. (Because updating a reference needs a GC write barrier due to the generational garbage collector.)

So, you're comparing two functions that differ both in algorithm an implementation technique. You should rewrite the second one not to use references, that way you can compare just algorithms. Or rewrite the first one to use references. But why use references when it's the wrong thing? :)

You can compare the algorithmic complexities.

The first is O(1);

the second is O(n)