fibonacci memoization speed tests

solrize 2012-06-24 23:39:56.134007 UTC

1{- benchmark of a bunch of memoization strategies for fibonacci
2 "Fibonacci is the E. Coli of Haskell optimization experiments"
3 (misquoting someone)
4
5 Author: solrize, June 24, 2012
6
7 Inspiration:
8 http://stackoverflow.com/questions/4980102/how-does-data-memocombinators-work
9
10 Output on 2.5 ghz Intel Core 2, GHC 7.03 with -O2:
11 name usec (N = 50000)
12 MT 484926
13 linear 292955
14 array 180973
15 conal 534918
16 luqui 939857
17 direct 94985
18 intmap 371944
19 ugly 528920
20 STArray 327949
21-}
22{-# LANGUAGE BangPatterns, RankNTypes #-}
23
24module Main (main) where
25
26import System.CPUTime
27import Data.Array
28import qualified Data.MemoTrie as Conal
29import qualified Data.MemoCombinators as Luqui
30import GHC.IO
31import Text.Printf
32import qualified Data.IntMap as IM
33import Control.Monad.ST
34import Data.STRef
35import qualified Data.MemoUgly as Ugly
36import Data.Array.ST
37import Data.Maybe (fromJust)
38
39testValue = 50000
40
41type Fibn = Int -> Integer
42
43data MTree a r = MTree {
44 v0 :: r
45 , v1 :: r
46 , m0 :: MTree a r
47 , m1 :: MTree a r
48 } deriving Show
49
50memo :: Integral a => (a -> r) -> (a -> r)
51memo f = tfind (r 1 0)
52 where
53 r p k = MTree { v0 = f k, v1 = f (k+p)
54 , m0 = r (p*2) k, m1 = r (p*2) (k+p)
55 }
56
57 tfind (MTree {v0=v0, v1=v1, m0=m0, m1=m1}) k
58 | k < 0 = error "k<0 is left as an exercise"
59 | k == 0 = v0
60 | k == 1 = v1
61 | even k = tfind m0 k2
62 | odd k = tfind m1 k2
63 where
64 k2 = k`div`2
65
66-- direct tail recursion
67fibDirect :: Fibn
68fibDirect n = go n 0 1
69 where
70 go 0 a b = a
71 go !n !a !b = go (n-1) b (a+b)
72
73-- MTree infinite tree memoization
74fibMT :: Fibn
75fibMT = memo fib0
76 where
77 fib0 0 = 0
78 fib0 1 = 1
79 fib0 n = fibMT (n-1) + fibMT (n-2)
80
81-- reference linear implementation
82fibLinear :: Fibn
83fibLinear n = fibLinear' 0 1 !! n
84 where
85 fibLinear' m n = m : fibLinear' n (m+n)
86
87-- array-based memoization
88memarray = array (0,testValue) [(i,fibArray i)|i<-[0..testValue]]
89fibArray :: Fibn
90fibArray 0 = 0
91fibArray 1 = 1
92fibArray n = memarray ! (n-1) + memarray ! (n-2)
93
94-- Conal Eliot's MemoTrie
95fibConal :: Fibn
96fibConal = Conal.memo fibr
97 where
98 fibr 0 = 0
99 fibr 1 = 1
100 fibr n = fibConal (n-1) + fibConal (n-2)
101
102-- Luqui's memo combinators
103fibLuqui :: Fibn
104fibLuqui = Luqui.integral fibr
105 where
106 fibr 0 = 0
107 fibr 1 = 1
108 fibr n = fibLuqui (n-1) + fibLuqui (n-2)
109
110-- imperative-ish manual intmap memoization
111fibIntMap :: Fibn
112fibIntMap n = runST $ do
113 xmemo <- newSTRef $ IM.fromList [(0,0),(1,1)]
114 let fibx n = do
115 v <- readSTRef xmemo
116 case IM.lookup n v of
117 Just x -> return x
118 Nothing -> do
119 a <- fibx (n-1)
120 b <- fibx (n-2)
121 let c = a + b
122 modifySTRef xmemo (IM.insert n c)
123 return c
124 fibx n
125
126-- Uglymemo (uses unsafePerformIO)
127fibUgly :: Fibn
128fibUgly = Ugly.memo fibr
129 where
130 fibr 0 = 0
131 fibr 1 = 1
132 fibr n = fibUgly (n-1) + fibUgly (n-2)
133
134-- mutable array per jmcarthur's suggestion of a hash table
135fibSTArray :: Fibn
136fibSTArray n = imemo
137 where
138 imemo = runST $ do
139 xmemo <- newArray (0, testValue) Nothing
140 :: ST s (STArray s Int (Maybe Integer))
141 writeArray xmemo 0 (Just 0)
142 writeArray xmemo 1 (Just 1)
143 let fibr k = do
144 r <- readArray xmemo k
145 case r of
146 Just x -> return x
147 Nothing -> do
148 a <- fibr (k-1)
149 b <- fibr (k-2)
150 let c = a+b
151 writeArray xmemo k (Just c)
152 return c
153 fibr n
154{-
155test = a1 == a2 && a2 == a3 where
156 a1 = fibMT k
157 a2 = fibLinear k
158 a3 = fibArray k
159 k = 12345
160-}
161
162
163main = do
164 putStrLn $ "name usec (N = " ++ show testValue ++ ")"
165 p "MT" fibMT
166 p "linear" fibLinear
167 p "array" fibArray
168 p "conal" fibConal
169 p "luqui" fibLuqui
170 p "direct" fibDirect
171 p "intmap" fibIntMap
172 p "ugly" fibUgly
173 p "STArray" fibSTArray
174
175 where
176 p label f = do
177 t0 <- getCPUTime
178 a <- evaluate (f testValue)
179 t1 <- getCPUTime
180 let dt = (t1-t0)`div`(10^6) -- microseconds
181 printf "%-8s %7d\n" label dt