speeding up fibonacci with memoizing
I just thought this was interesting, so I would share it. Thanks to whoever it was on #haskell who helped me, sorry I can't remember who. -- horribly slow. try slow_fibs 30, not too much higher than that and it hangs slow_fibs = map slow_fib [1..] slow_fib 1 = 1 slow_fib 2 = 1 slow_fib n = ( slow_fib (n-2) ) + (slow_fib(n-1)) -- versus, try memoized_fibs !! 10000 memoized_fibs = map memoized_fib [1..] memoized_fib = ((map fib' [0 ..]) !!) where fib' 0 = 0 fib' 1 = 1 fib' n = memoized_fib (n - 1) + memoized_fib (n - 2)
On 2/18/07, Thomas Hartman
-- versus, try memoized_fibs !! 10000 memoized_fibs = map memoized_fib [1..] memoized_fib = ((map fib' [0 ..]) !!) where fib' 0 = 0 fib' 1 = 1 fib' n = memoized_fib (n - 1) + memoized_fib (n - 2)
How about this: zip_fibs = 0 : 1 : [x+y | ~(x,y) <- zip zip_fibs (tail zip_fibs)] zip_fib = (!!) zip_fibs A naïve performance test: ---------------------------------------------------------------- module Main (main) where import System (getArgs) zip_fibs = 0 : 1 : [x+y | ~(x,y) <- zip zip_fibs (tail zip_fibs)] zip_fib = (!!) zip_fibs memoized_fibs = map memoized_fib [1..] memoized_fib = ((map fib' [0 ..]) !!) where fib' 0 = 0 fib' 1 = 1 fib' n = memoized_fib (n - 1) + memoized_fib (n - 2) main = do [func,arg] <- getArgs let n = read arg case func of "zip_fib" -> do putStrLn $ "Using 'zip_fib' to calculate fib of " ++ arg putStrLn $ " -> " ++ show (zip_fib n) "memoized_fib" -> do putStrLn $ "Using 'memoized_fib' to calculate fib of " ++ arg putStrLn $ " -> " ++ show (memoized_fibs !! n) ---------------------------------------------------------------- $ rm *.o $ ghc --version The Glorious Glasgow Haskell Compilation System, version 6.4.2 $ ghc -O fib.hs $ time ./a.out zip_fib 20000 Using 'zip_fib' to calculate fib of 20000 -> 25311623237323612422[...]39857683971213093125 real 0m0.140s user 0m0.106s sys 0m0.001s $ time ./a.out memoized_fib 20000 Using 'memoized_fib' to calculate fib of 20000 -> 25311623237323612422[...]39857683971213093125 real 0m30.787s user 0m29.509s sys 0m0.156s $ time ./a.out zip_fib 200000 Using 'zip_fib' to calculate fib of 200000 -> 15085683557988938992[...]52246259408175853125 real 0m24.790s user 0m23.649s sys 0m0.159s No, I *won't* try './a.out memoized_fib 200000' ;-). -- Felipe.
Besides memoizing, you might want to use the fact that: fib (2*k) == (fib (k+1))^2 - (fib (k-1))^2 fib (2*k-1) == (fib k)^2 + (fib (k-1))^2 Regards, Yitz
On 2/18/07, Yitzchak Gale
Besides memoizing, you might want to use the fact that:
fib (2*k) == (fib (k+1))^2 - (fib (k-1))^2 fib (2*k-1) == (fib k)^2 + (fib (k-1))^2
Nice one: $ time ./a.out another_fib 200000 Using 'another_fib' to calculate fib of 200000 -> 15085683557988938992[...]52246259408175853125 real 0m1.177s user 0m1.051s sys 0m0.017s Implementation details: ----------------------------------------- another_fibs = 0 : 1 : 1 : map f [3..] where square x = x * x sqfib = square . another_fib 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 another_fib = (!!) another_fibs ----------------------------------------- Conclusion: it's often better to improve the algorithm than the implementation =). -- Felipe.
G'day all.
On 2/18/07, Yitzchak Gale
Besides memoizing, you might want to use the fact that:
fib (2*k) == (fib (k+1))^2 - (fib (k-1))^2 fib (2*k-1) == (fib k)^2 + (fib (k-1))^2
Quoting Felipe Almeida Lessa
Implementation details: ----------------------------------------- another_fibs = 0 : 1 : 1 : map f [3..] where square x = x * x sqfib = square . another_fib 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 another_fib = (!!) another_fibs -----------------------------------------
First off, your memo structure is a linked list, which takes O(n) time to access the nth element. The first call takes O(n) time, the second takes O(n/2) time, the third takes O(n/4) time etc etc, so in total, it's O(n). That's the same complexity as the naive iterative algorithm: fib n = fib' 0 1 !! n where fib' a b = a : fib' b (a+b) Secondly, the memo data structure here leaks memory. If you need fib 2000000 once and only low Fibonacci numbers after that, you keep the data structure up to size 2000000, even though you don't need it. Now that's fine for a benchmark, but you should never do this in a library. Taking one and two together, it seems that it would be better to use an array. Let's try that: fib :: Integer -> Integer fib n | n < 0 = error "fib" -- A slight lie, fixed below. | n < fromIntegral memoSize = memoTable ! fromIntegral n | even n = let n2 = n `div` 2 a = fib (n2+1) b = fib (n2-1) in a*a - b*b | otherwise = let n2 = (n+1) `div` 2 a = fib n2 b = fib (n2-1) in a*a + b*b where memoSize :: Int memoSize = 10000 memoTable = array (0,memoSize-1) (take memoSize (fibs 0 1)) where fibs a b = a : fibs b (a+b) That keepe the memory leak under control, but there's another problem. If n >= memoSize, then there will be two recursive calls spawned, which will spawn two recursive calls... The complexity is O(2^n), which is exactly the same problem as the naive recursive Fibonacci implementation. We've effectively just made the constant factor extremely low by optimising a huge number of base cases. Now consider the recursive calls. The even case needs fib (n2-1) and fib (n2+1), and the odd case needs fib n2 and fib (n2-1). If we have a = fib n2 and b = fib (n2-1), we can trivially compute fib (n2+1); it's just a+b. So let's just modify the recursive call to return two adjacent Fibonacci numbers. That way, we only have one recursive call, and the overall complexity should be O(log n). We want something like this: fib :: Integer -> Integer fib n | n < 0 -- Because fib (n+1) = fib n + fib (n-1), we can extend -- Fibonacci numbers below 0. = let n' = -n in if even n then -fib n' else fib n' | otherwise = fst (fib' n) fib' :: Integer -> (Integer,Integer) The base cases: fib' n | n < fromIntegral memoSize = memoTable ! fromIntegral n where memoSize :: Int memoSize = 10000 memoTable = listArray (0,memoSize-1) (take memoSize (fibs 0 1)) where fibs a b = (a,b) : fibs b (a+b) The array contains the pairs (0,1), (1,1), (1,2), (2,3), (3,5) etc. If you think about it, this is redundant; we're holding twice the number of Integers that we need to. So let's optimise that a bit: fib' n | q < fromIntegral memoSize = case memoTable ! fromIntegral q of p@(a,b) | r == 0 -> p | otherwise -> (b, a+b) 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) Finally, the recursive case. A little arithmetic gives: fib' n | q < fromIntegral memoSize = {- as before -} | r == 0 = let (a,b) = fib' (q-1) c = a+b c2 = c*c in (c2 - a*a, c2 + b*b) | otherwise = let (a,b) = fib' q c = a+b a2 = a*a in (b*b + a2, c*c - a2) where {- as before -} Now that's an industrial-strength Fibonacci. It's O(log n) not including the cost of adding and multiplying large Integers, and uses a bounded amount of memory between calls, making it suitable for a library. The slowest part of the test program is actually the bit that prints the number. So I used this driver program: main :: IO () main = do (n:_) <- getArgs putStrLn "another_fib" putStrLn (another_fib (read n) `seq` "done") I also used a slow machine so we can see the difference. Your version: % time ./fibtest 200000 another_fib done real 0m1.173s user 0m1.041s sys 0m0.121s And mine: % time ./fibtest 200000 fib done real 0m0.312s user 0m0.289s sys 0m0.021s
Conclusion: it's often better to improve the algorithm than the implementation =).
And it's even better to do both. Cheers, Andrew Bromage
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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@stefans:/tmp$
Would someone please update the entries on our 'archive of fibs' page? http://www.haskell.org/haskellwiki/The_Fibonacci_sequence Cheers. stefanor:
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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@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@stefans:/tmp$ _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
G'day all.
Quoting Stefan O'Rear
Prior art trumps all. (by a few %) granted it doesn't do much memoizing anymore :)
Ah, butbutbut... of course the Gosper/Salamin one is going to be faster if you only compute one Fibonacci number per instance. The memoed version is optimised for programs that want more than one. Cheers, Andrew Bromage
Stefan O'Rear wrote:
Prior art trumps all. (by a few %) granted it doesn't do much memoizing anymore :)
gs > ajb > f > d > u, it, z > s > n
[snip] Nice. I took the opportunity to polish my generic linear recurrence module somewhat and test its speed. It does quite well I think: using http://int-e.home.tlink.de/haskell/LinRec.hs and defining
import qualified LinRec as L
-- generic linear recurrence generator genfib = L.get [1,1] [0,1]
I get: # ./t.sh gen 50000000 11.65user 0.06system 0:11.71elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k 0inputs+0outputs (0major+17092minor)pagefaults 0swaps # ./t.sh gs 50000000 4.67user 0.06system 0:04.79elapsed 98%CPU (0avgtext+0avgdata 0maxresident)k 0inputs+0outputs (0major+11746minor)pagefaults 0swaps for a slowdown of about 2.5. Bertram
On Sun, Feb 18, 2007 at 06:15:25PM -0500, ajb@spamcop.net wrote:
Now that's an industrial-strength Fibonacci. It's O(log n) not including the cost of adding and multiplying large Integers, and uses a bounded amount of memory between calls, making it suitable for a library. The slowest part of the test program is actually the bit that prints the number. So I used this driver program:
I've been here before. http://www.haskell.org/pipermail/haskell-cafe/2005-January/008839.html -- wli
On Sun, 18 Feb 2007, Yitzchak Gale wrote:
Besides memoizing, you might want to use the fact that:
fib (2*k) == (fib (k+1))^2 - (fib (k-1))^2 fib (2*k-1) == (fib k)^2 + (fib (k-1))^2
Or, you know, go straight to the closed form for the fibonacci numbers! :) -- Mikael Johansson | To see the world in a grain of sand mikael@johanssons.org | And heaven in a wild flower http://www.mikael.johanssons.org | To hold infinity in the palm of your hand | And eternity for an hour
On Mon, Feb 19, 2007 at 08:47:39AM +0100, Mikael Johansson wrote:
On Sun, 18 Feb 2007, Yitzchak Gale wrote:
Besides memoizing, you might want to use the fact that:
fib (2*k) == (fib (k+1))^2 - (fib (k-1))^2 fib (2*k-1) == (fib k)^2 + (fib (k-1))^2
Or, you know, go straight to the closed form for the fibonacci numbers! :)
That's fine in the blessed realm of arithmatic rewrite rules, but here we need bitstrings, and computing large powers of irrational numbers is not exactly fast. Phi is definable in finite fields (modular exponentiation yay!) but modular-ation seems ... problematic. I have a gut feeling the p-adic rationals might help, but insufficient knowledge to formulate code. The GMP fibbonacci implementation is of my quasilinear recurrence family, not closed form. And lest we forget the obvious - by far the fastest way to implement fib in GHC Haskell: {-# OPTIONS_GHC -O2 -cpp -fglasgow-exts #-} module #ifdef fibimpl Main(main) #else Fibs #endif where import System.Environment import Array import List(unfoldr) #ifdef __GLASGOW_HASKELL__ import System.IO.Unsafe import Foreign.Marshal.Alloc import Foreign.Storable import GHC.Exts import Foreign.C.Types #endif -- same as before #ifdef __GLASGOW_HASKELL__ foreign import ccall "gmp.h __gmpz_fib_ui" _gfib :: Ptr Int -> CULong -> IO () foreign import ccall "gmp.h __gmpz_init" _ginit :: Ptr Int -> IO () gmpfib :: Int -> Integer gmpfib n = unsafePerformIO $ allocaBytes 12 $ \p -> do _ginit p _gfib p (fromIntegral n) I# sz <- peekElemOff p 1 I# pt <- peekElemOff p 2 return (J# sz (unsafeCoerce# (pt -# 8#))) #endif -- same as before stefan@stefans:~/fibbench$ ./h gs 100000000 12.84user 0.24system 0:13.08elapsed 100%CPU (0avgtext+0avgdata 0maxresident)k 0inputs+0outputs (0major+21082minor)pagefaults 0swaps stefan@stefans:~/fibbench$ ./h gmp 100000000 9.12user 0.42system 0:09.58elapsed 99%CPU (0avgtext+0avgdata 0maxresident)k 0inputs+0outputs (0major+35855minor)pagefaults 0swaps stefan@stefans:~/fibbench$
"Thomas Hartman"
module MemoFib where
The unexciting version
naive_fib 0 = 1 naive_fib 1 = 1 naive_fib n = naive_fib (n-1) + naive_fib (n-2)
The memoised version using a memoising fixpoint operator
fibonacci = memoFix fib where fib fib 0 = 1 fib fib 1 = 1 fib fib n = fib (n-1) + fib (n-2)
I suppose if you want to "put it in a library", you should just put fib in, and allow the user to call memoFix fib to make a new version when necessary? A memoising fixpoint operator. It works by putting the result of the first call of the function for each natural number into a data structure and using that value for subsequent calls ;-)
memoFix f = mf where memo = fmap (f mf) (naturals 1 0) mf = (memo !!!)
A data structure with a node corresponding to each natural number to use as a memo.
data NaturalTree a = Node a (NaturalTree a) (NaturalTree a)
Map the nodes to the naturals in this order: 0 1 2 3 5 4 6 7 ... Look up the node for a particular number
Node a tl tr !!! 0 = a Node a tl tr !!! n | odd n = tl !!! top | otherwise = tr !!! (top-1) where top = n `div` 2
We surely want to ba able to map on these things...
instance Functor NaturalTree where fmap f (Node a tl tr) = Node (f a) (fmap f tl) (fmap f tr)
If only so that we can write cute, but inefficient things like the below, which is just a NaturalTree such that naturals!!!n == n: naturals = Node 0 (fmap ((+1).(*2)) naturals) (fmap ((*2).(+1)) naturals) The following is probably more efficient (and, having arguments won't hang around at top level, I think) -- have I put more $!s than necessary?
naturals r n = Node n ((naturals $! r2) $! (n+r)) ((naturals $! r2) $! (n+r2)) where r2 = 2*r
Of course, if you want to take advantage of the pseudo O(n) lookup time of arrays, you could use a NaturalTree of arrays of some fixed size -- but arrays are O(log n) really... -- Jón Fairbairn Jon.Fairbairn@cl.cam.ac.uk
fib :: Integer -> Integer fib n = fibaux n 0 1 1 where fibaux :: Integer -> Integer -> Integer -> Integer -> Integer fibaux i a b c | i==0 = a | i/=0 = fibaux (i-1) b c (b+c)
Throwing in a trace statement in fibaux tells me that fibaux i a b c is not being memoized. If I do map fib [7..9], fibaux counts down to 0 afresh for each of 7, 8, and 9. Ideally, in map fib [0..n], fib (i-2) and fib (i-1) should be memoized and fib i would be evaluated in constant time. This is what happens if the loop is unrolled explicitly. marnes wrote:
fib :: Integer -> Integer fib n = fibaux n 0 1 1 where fibaux :: Integer -> Integer -> Integer -> Integer -> Integer fibaux i a b c | i==0 = a | i/=0 = fibaux (i-1) b c (b+c)
_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
If people write any new variants, please add them to: http://haskell.org/haskellwiki/The_Fibonacci_sequence :) westondan:
Throwing in a trace statement in fibaux tells me that fibaux i a b c is not being memoized. If I do map fib [7..9], fibaux counts down to 0 afresh for each of 7, 8, and 9. Ideally, in map fib [0..n], fib (i-2) and fib (i-1) should be memoized and fib i would be evaluated in constant time. This is what happens if the loop is unrolled explicitly.
marnes wrote:
fib :: Integer -> Integer fib n = fibaux n 0 1 1 where fibaux :: Integer -> Integer -> Integer -> Integer -> Integer fibaux i a b c | i==0 = a | i/=0 = fibaux (i-1) b c (b+c)
_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
On Tue, 6 Nov 2007, marnes wrote:
fib :: Integer -> Integer fib n = fibaux n 0 1 1 where fibaux :: Integer -> Integer -> Integer -> Integer -> Integer fibaux i a b c | i==0 = a | i/=0 = fibaux (i-1) b c (b+c)
participants (14)
-
ajb＠spamcop.net -
Bertram Felgenhauer -
Dan Weston -
Don Stewart -
dons＠cse.unsw.edu.au -
Felipe Almeida Lessa -
Henning Thielemann -
Jón Fairbairn -
marnes -
Mikael Johansson -
Stefan O'Rear -
Thomas Hartman -
William Lee Irwin III -
Yitzchak Gale