[Haskell-cafe] speeding up fibonacci with memoizing

Sun Feb 18 22:45:09 EST 2007

Prior art trumps all.  (by a few %)  granted it doesn't do much memoizing anymore :)

gs > ajb > f > d > u, it, z > s > n

stefan at stefans:/tmp$ ./h n 42
28.92user 0.14system 0:29.85elapsed 97%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+494minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h d 42
0.00user 0.00system 0:00.00elapsed 0%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+254minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h z 100
0.00user 0.00system 0:00.00elapsed 200%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+259minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h z 10000
0.03user 0.00system 0:00.03elapsed 105%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+746minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h z 100000
3.46user 0.02system 0:03.48elapsed 100%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+1981minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h d 100000
1.00user 0.00system 0:01.01elapsed 100%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+759minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h s 100000
3.70user 0.03system 0:03.73elapsed 100%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+2175minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h it 100000
3.43user 0.02system 0:03.46elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+1981minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h u 100000
3.41user 0.03system 0:03.45elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+1981minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h u 200000
17.34user 0.05system 0:17.44elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+3200minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h it 200000
17.38user 0.06system 0:18.99elapsed 91%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+3199minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h it 200000
17.31user 0.06system 0:17.70elapsed 98%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+3200minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h z 200000
17.34user 0.07system 0:17.42elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+3199minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h s 200000
20.15user 0.09system 0:20.25elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+3591minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h d 200000
4.20user 0.02system 0:04.24elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+758minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h d 100000
1.02user 0.01system 0:01.03elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+758minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h f 200000
0.12user 0.02system 0:00.14elapsed 102%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+2301minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h f 1000000
0.64user 0.08system 0:00.72elapsed 100%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+8456minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h f 3000000
2.58user 0.38system 0:02.96elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+33037minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h f 5000000
3.46user 0.40system 0:03.87elapsed 100%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+33036minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h ajb 5000000
0.52user 0.02system 0:00.54elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+2181minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h gs 5000000
0.39user 0.01system 0:00.41elapsed 100%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+1747minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h gs 50000000
5.85user 0.11system 0:05.96elapsed 100%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+11183minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h f 10000000
6.93user 0.91system 0:07.95elapsed 98%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+66059minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h gs 10000000
0.90user 0.04system 0:00.97elapsed 98%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+3379minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h ajb 10000000
1.08user 0.04system 0:01.12elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+3584minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h ajb 100000000
14.09user 0.25system 0:14.42elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+23586minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h gs 100000000
13.17user 0.20system 0:13.48elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+19588minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h ajb 300000000
49.05user 0.80system 0:50.71elapsed 98%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+64948minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h gs 300000000
46.21user 0.60system 0:46.89elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+55454minor)pagefaults 0swaps
stefan at stefans:/tmp$ ./h gs 300000000
46.29user 0.67system 0:47.62elapsed 98%CPU (0avgtext+0avgdata 0maxresident)k
0inputs+0outputs (0major+57529minor)pagefaults 0swaps
stefan at stefans:/tmp$ cat h
#! /bin/sh
ghc-6.6 -v0 -fforce-recomp -O2 -cpp -Dfibimpl=${1}fib --make Y.hs
/usr/bin/time ./Y $2
stefan at stefans:/tmp$ cat Y.hs
{-# OPTIONS_GHC -O2 -cpp #-}
module
#ifdef fibimpl
Main(main)
#else
Fibs
#endif
where
import System
import Array
import List(unfoldr)

fix f = let x = f x in x -- is this in h98?

--- Exponential fibs

-- *The* naive fib
nfib :: Int -> Integer
nfib 0 = 0
nfib 1 = 1
nfib n = nfib (n-1) + nfib (n-2)

--- Quadratic fibs

-- Zipping fib
zfib = (!!) zfibs where zfibs :: [Integer]
                        zfibs = 0 : 1 : zipWith (+) zfibs (tail zfibs)

-- Scanning fib
sfib = (!!) (fix ((0:) . scanl (+) (1 :: Integer)))

-- Iterative list fib, explicit recursion
itfib = (!!) (ifibs 0 1)
    where ifibs :: Integer -> Integer -> [Integer]
          ifibs a b = a : ifibs b (a+b)

-- Iterative list fib, unfoldr
ufib = (!!) (unfoldr fibit (0,1))
    where fibit :: (Integer,Integer) -> Maybe (Integer, (Integer,Integer))
          fibit (a,b) = Just (a, (b, a+b))

-- Iterative fib, explicitly deforested
dfib n = dfibl n 0 1 where dfibl :: Int -> Integer -> Integer -> Integer
                           dfibl 0 a b = b `seq` a
                           dfibl n a b = n `seq` a `seq` b `seq` dfibl (n-1) b (a+b)

--- Quasilinear fibs

-- Gosper/Salamin fib, as seen in HAKMEM #12
gsfib' :: Int -> (Integer,Integer)
gsfib' n | n == 0 = (0,1)
         | odd  n = case gsfib' (n-1)       of (a,b) -> (a+b, a)
         | even n = case gsfib' (n `div` 2) of (a,b) -> (a*(a+b+b), a*a + b*b)
gsfib n = fst (gsfib' n)

-- Felipe's fib, using a recurrence relation + memoization
ffibs :: [Integer]
ffibs = 0 : 1 : 1 : map f [3..]
     where
       square x = x * x
       sqfib = square . ffib
       f n | even n = sqfib (k+1) - sqfib (k-1) where k = n `div` 2
       f n          = sqfib k + sqfib (k-1) where k = (n + 1) `div` 2
ffib = (!!) ffibs

-- AJB's super-optimized memoizing fib
ajbfib = fst . ajbfib'
ajbfib' :: Int -> (Integer,Integer)
ajbfib' n
    | q < fromIntegral memoSize
      = case memoTable ! fromIntegral q of
          p@(a,b) | r == 0    -> p
                  | otherwise -> (b, a+b)
    | r == 0
      = let (a,b) = ajbfib' (q-1)
            c = a+b
            c2 = c*c
        in (c2 - a*a, c2 + b*b)
    | otherwise
      = let (a,b) = ajbfib' q
            c = a+b
            a2 = a*a
        in (b*b + a2, c*c - a2)
    where
        (q,r) = n `divMod` 2

        memoSize :: Int
        memoSize = 10000

        memoTable
          = listArray (0,memoSize-1) (take memoSize (fibs 0 1))
          where
            fibs a b = (a,b) : let ab = a+b in fibs ab (ab+b)

#ifdef fibimpl
main = do [ x ] <- mapM readIO =<< System.getArgs
          seq (fibimpl x) (return ())
#endif
stefan at stefans:/tmp$