👨‍💻 about me home CV/Resume 🖊️ Contact Github LinkedIn I’m a Haskeller 📝 Blog Freedom, privacy, tutorials… 🏆 Best of Fizzbuzz makex LuaX Calculadoira panda upp Haskell todo pwd TPG Nextcloud Git

🆕 since August 2022: LuaX is a Lua eXtended interpretor/cross compiler providing a bunch of useful modules (statically linked, no dependency). Nice integration with upp (new functions and modules available to extend upp macros) and also with a soon released but yet confidential project about actor oriented programming!
💣 Kick GAFAMs out (✔️ ǝlƃooפ, ✔️ ʞooqǝɔɐℲ, ✔️ uozɐɯ∀): Stop giving our soul and money to evils, be free and respectful!
📰 Friday 2. April 2021: upp is a panda companion. It’s a Lua-scriptable lightweight text preprocessor.
🆕 since December 2020: Playing with the actor model in an embedded multicore context. C imperative components become C stream pure functions with no side effect ➡️ C low level programming with high level pure functional programming properties 🏆
📰 Saturday 30. January 2021: Playing with Pandoc Lua filters in Lua. panda is a lightweight alternative to abp providing a consistent set of Pandoc filters (text substitution, file inclusion, diagrams, scripts, …).
🆕 Sunday 24. May 2020: Working at EasyMile for more than 5 years. Critical real-time software in C, simulation and monitoring in Haskell ➡️ perfect combo! It’s efficient and funny ;-)

Super fast recursive Fibonacci implementation

Christophe Delord

11 June 2022

License

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

Introduction

The Fibonacci sequence is a classical hello-world application for functional programming. The challenge here is to get a fast implementation. We will first show two classical implementations: the trivial recursive definition that is very slow and the iterative version that is slightly faster.

And a third one that is blazing fast.

Definition

The Fibonacci sequence is defined recursively:

\[ \begin{align} F_0 &= 1 \\ F_1 &= 1 \\ F_n &= F_{n-1} + F_{n-2} \end{align} \]

\(n\) \(F_n\) \(F_n = F_{n-1} + F_{n-2}\)
0 1
1 1
2 2 \(F_{2} = F_{1} + F_{0}\)
3 3 \(F_{3} = F_{2} + F_{1}\)
4 5 \(F_{4} = F_{3} + F_{2}\)
5 8 \(F_{5} = F_{4} + F_{3}\)
6 13 \(F_{6} = F_{5} + F_{4}\)
7 21 \(F_{7} = F_{6} + F_{5}\)
8 34 \(F_{8} = F_{7} + F_{6}\)
9 55 \(F_{9} = F_{8} + F_{7}\)
10 89 \(F_{10} = F_{9} + F_{8}\)

Trivial recursive Fibonacci sequence definition

A trivial implementation of this sequence in Haskell is:

recursiveFib :: Int -> Integer
recursiveFib 0 = 1
recursiveFib 1 = 1
recursiveFib n = recursiveFib (n-1) + recursiveFib (n-2)

The objective here is to compare implementations for large values of \(n\). For such values \(F_n\) can not be coded by machine integers. \(F_n\) can be an arbitrary large integer (Integer type in Haskell).

This function is very simple but also very slow. Its complexity is \(\mathcal{O}(F_n)\) and \(F_n\) increases very fast since \(F_n \sim \frac{1}{\sqrt{5}}\left(\frac{1+\sqrt{5}}{2}\right)^n\).

Iterative Fibonacci definition

\(F_n\) can also be constructed iteratively from \(F_0\) to \(F_n\):

\[ \begin{align} F_0 &= 1 \\ F_1 &= 1 \\ F_2 &= F_1 + F_0 \\ F_3 &= F_2 + F_1 \\ ... \\ F_n &= F_{n-1} + F_{n-2} \\ \end{align} \]

iterativeFib :: Int -> Integer
iterativeFib n = loop n 1 1
    where
        loop 0 a b = b
        loop n a b = loop (n-1) (a+b) a

This function is faster since it avoids computing Fibonacci terms several times. Its complexity is \(\mathcal{O}(n)\).

Fast Fibonacci implementation

We can exploit some properties of the Fibonacci sequence to improve its computation.

Let’s try to apply the definition several times to compute several steps at once:

\[ \begin{align} F_n &= F_{n-1} &+ F_{n-2} \\ &= F_{n-2} + F_{n-3} &+ F_{n-2} \\ &= 2 F_{n-2} &+ F_{n-3} \\ &= 2 \left( F_{n-3} + F_{n-4} \right) &+ F_{n-3} \\ &= 3 F_{n-3} &+ 2 F_{n-4} \\ &= 3 \left( F_{n-4} + F_{n-5} \right) &+ 2 F_{n-4} \\ &= 5 F_{n-4} &+ 3 F_{n-5} \\ \end{align} \]

The coefficients seem to be Fibonacci numbers.

Let’s prove by induction (over \(k\)) that:

\[ F_n = F_k F_{n-k} + F_{k-1} F_{n-k-1} \]

The property is obviously true for \(k = 1\):

\[ \begin{align} F_n &= F_1 F_{n-1} &+ F_{1-1} F_{n-1-1} \\ &= F_1 F_{n-1} &+ F_0 F_{n-2 } \\ &= 1 F_{n-1} &+ 1 F_{n-2 } \\ &= F_{n-1} &+ F_{n-2 } \\ \end{align} \]

Lets prove that \(F_n = F_k F_{n-k} + F_{k-1} F_{n-k-1} \implies F_n = F_{k+1} F_{n-k-1} + F_k F_{n-k-2}\).

\[ \begin{align} F_{k+1} F_{n-k-1} + F_k F_{n-k-2} &= \left( F_k + F_{k-1} \right) F_{n-k-1} + F_k F_{n-k-2} \\ &= F_k \left( F_{n-k-1} + F_{n-k-2} \right) + F_{k-1} F_{n-k-1} \\ &= F_k F_{n-k} + F_{k-1} F_{n-k-1} \\ \end{align} \]

The Fibonacci sequence can then be defined by:

\[ \begin{align} F_0 &= 1 \\ F_1 &= 1 \\ F_n &= F_k F_{n-k} + F_{k-1} F_{n-k-1}, \forall n \ge 2, \forall k \in [1, n-1] \end{align} \]

The trivial implementation is slow because of the double recursivity. We now have four recursions but:

  1. The call tree is reduced by jumping from \(n\) to \(n-k\) instead of just \(n-1\)
  2. By choosing \(k\) wisely some redundant computations can be avoided

Intuitively if \(k\) is close to \(n-k\) then subtrees will be close too. If \(k = n-k\) (i.e. \(k = \frac{n}{2}\)) then the first term of \(F_n\) is a square (\(F_k = F_{n-k}\)).

First case: n is even (\(n = 2k\))

\[ \begin{align} F_n &= F_k F_{n-k} &+ F_{k-1} F_{n-k-1} \\ &= F_k F_k &+ F_{k-1} F_{k-1} &&(n-k-1 = 2k-k-1 = k-1) \\ &= F_k^2 &+ F_{k-1}^2 \\ \end{align} \]

Second case: n is odd (\(n = 2k+1\))

\[ \begin{align} F_n &= F_k F_{n-k} &+ F_{k-1} F_{n-k-1} \\ &= F_k F_{k+1} &+ F_{k-1} F_{n-k-1} && (n-k = 2k+1-k = k+1) \\ &= F_k F_{k+1} &+ F_{k-1} F_{k} && (n-k-1 = k) \\ &= F_k \left( F_k + F_{k-1} \right) &+ F_{k-1} F_{k} \\ &= F_k^2 + 2 F_k F_{k-1} \\ &= F_k \left(F_k + 2 F_{k-1} \right) \\ \end{align} \]

This definition can be easily implemented in Haskell:

fastFib :: Int -> Integer
fastFib 0 = 1
fastFib 1 = 1
fastFib n =
    case r of
        0 -> fk^2 + fk1^2
        1 -> fk * (fk + 2*fk1)
    where
        (k, r) = n `divMod` 2
        fk = fastFib k
        fk1 = fastFib (pred k)

This function is way faster than the previous ones and can deal with very large numbers. Its complexity is \(\mathcal{O}(F_{\log_2 n})\) (I have no proof yet, just an intuition because the depth of the tree is divided by two every steps and there is still a double recursion).

Tests

Unit tests

The three functions shall return the same values:

fibs = [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597]

testFib :: String -> (Int -> Integer) -> IO ()
testFib name f = do
    forM_ (zip3 [0..] fibs (map f [0..])) $ \(i, ref, fi) ->
        when (fi /= ref) $
            error ("ERROR: "++name++" "++show i++" = "++show fi++" (should be "++show ref++")")
    putStrLn $ "    - "++name++": OK"

unitTests = do
    putStrLn "* Unit tests:"
    putStrLn ""
    testFib "recursiveFib" recursiveFib
    testFib "iterativeFib" iterativeFib
    testFib "fastFib" fastFib
    putStrLn ""

Performance tests

Lets compare the performances of the three function:

perfTests = do
    putStrLn "* Performance tests:"
    putStrLn ""
    putStrLn "  n   | recursiveFib | iterativeFib | fastFib "
    putStrLn "  ---:|-------------:|-------------:|--------:"
    forM_ [20..40] $ \i -> do
        t1 <- fmtTime <$> profile recursiveFib i
        t2 <- fmtTime <$> profile iterativeFib i
        t3 <- fmtTime <$> profile fastFib i
        putStrLn $ show i++" | "++t1++" | "++t2++" | "++t3
    forM_ [20000, 40000 .. 300000] $ \i -> do
        let t1 = ""
        t2 <- fmtTime <$> profile iterativeFib i
        t3 <- fmtTime <$> profile fastFib i
        putStrLn $ show i++" | "++t1++" | "++t2++" | "++t3
    forM_ [2500000, 5000000 .. 40000000] $ \i -> do
        let t1 = ""
        let t2 = ""
        t3 <- fmtTime <$> profile fastFib i
        putStrLn $ show i++" | "++t1++" | "++t2++" | "++t3
    putStrLn ""

profile :: (Int -> Integer) -> Int -> IO (Double, Integer)
profile f n = timeItT $ do
    let fn = f n
    fn `deepseq` (return ())
    return fn

fmtTime :: (Double, Integer) -> String
fmtTime (t, f) | t < 1e-6 = show (round (t*1e9)) ++ " ns"
fmtTime (t, f) | t < 1e-3 = show (round (t*1e6)) ++ " µs"
fmtTime (t, f) | t < 1e-0 = show (round (t*1e3)) ++ " ms"
fmtTime (t, f)            = show (round t) ++ " s"

Large values

largeFibValue = do
    let n = 100*1000
    let fn = fastFib n
    let fn' = fmt fn
    let nbDigits = length (filter (/=' ') fn')
    putStr "* Large number example: "
    putStrLn $ "fastFib "++fmt (fromIntegral n)++" = "++fn'++" ("++fmt (fromIntegral nbDigits)++" digits)"

fmt :: Integer -> String
fmt = intercalate " " . reverse . map reverse . chunksOf 3 . reverse . show

Tests results

Tests made on a Intel(R) Core(TM) i7-9700 CPU @ 3.00GHz powered by Fedora Linux 37 and The Glorious Glasgow Haskell Compilation System, version 9.2.5.