👨‍💻 about me home CV/Resume 🖊️ Contact Github LinkedIn I’m a Haskeller 📝 Blog Freedom, privacy, tutorials… 🏆 Best of panda upp Haskell abp pp hCalc bl lapp todo pwd TPG Nextcloud Git BitTorrent

💣 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 ;-)
🚌 And we are recruiting! Contact if you are interested in Haskell or embedded softwares (or both).

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 36 and The Glorious Glasgow Haskell Compilation System, version 8.10.7.