Fibonacci Numbers

July 30, 2010

One of the first functions taught to programmers who are just learning about recursion is the function to compute the fibonacci numbers. The naive function takes exponential time, as each recursive call must compute the values of smaller fibonacci numbers, so programmers are next taught how to remove the recursion by explicitly storing state information, giving a linear-time iterative algorithm. The upshot is that programming students are left with the impression that recursion is bad and iteration is good.

It is actually possible to improve the performance with a logarithmic-time algorithm. Consider the matrices

m = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}, m^2 = \begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix}, m^3 = \begin{pmatrix} 3 & 2 \\ 2 & 1 \end{pmatrix}, m^4 = \begin{pmatrix} 5 & 3 \\ 3 & 2 \end{pmatrix}, m^5 = \begin{pmatrix} 8 & 5 \\ 5 & 3 \end{pmatrix}.

Each time the matrix is multiplied by itself, the number in the lower left-hand corner is the next fibonacci number; for instance, F4=3 (F0=0 is a special case). Of course, powering can be done using a binary square-and-multiply algorithm, as in the ipow and expm functions of the Standard Prelude, giving a logarithmic algorithm for computing the nth fibonacci number.

Your task is to write the three fibonacci functions — exponential, linear, and logarithmic — described above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

11 Responses to “Fibonacci Numbers”

  1. [...] Praxis – Fibonacci Numbers By Remco Niemeijer In today’s Programming Praxis we have to provide three different methods of calculating the ever-popular [...]

  2. Remco Niemeijer said

    My Haskell solution (see http://bonsaicode.wordpress.com/2010/07/30/programming-praxis-fibonacci-numbers/ for a version with comments):

    import Data.List
    
    fibexp :: Int -> Integer
    fibexp 0 = 0
    fibexp 1 = 1
    fibexp n = fibexp (n - 1) + fibexp (n - 2)
    
    fiblin :: Int -> Integer
    fiblin n = fibs !! n where fibs = 0:1:zipWith (+) fibs (tail fibs)
    
    mult :: Num a => [[a]] -> [[a]] -> [[a]]
    mult a b = [map (sum . zipWith (*) r) $ transpose b | r <- a]
    
    matrixpower :: [[Integer]] -> Int -> [[Integer]]
    matrixpower m 1 = m
    matrixpower m n = (if even n then id else mult m) $
                      matrixpower (mult m m) (div n 2)
    
    fiblog :: Int -> Integer
    fiblog 0 = 0
    fiblog n = matrixpower [[1,1],[1,0]] n !! 1 !! 0
    
  3. Graham said

    I cheated a bit to get the logic of my matrix_mul and matrix_pow right, looking at Remco’s very clean answer and some other sources (I haven’t had coffee yet today…)

    Here’s an inelegant Python implementation; apologies for the ugly matrix_mul definition.

    def fib_exp(n):
        if n < 2:
            return 1
        else:
            return fib_exp(n-1) + fib_exp(n-2)
    
    def fib_lin(n):
        c = n
        a, b = 1, 1
        while c > 0:
            a, b, c = b, a + b, c - 1
        return a
    
    def matrix_mul(a, b):
        return [[a[0][0]*b[0][0]+a[0][1]*b[1][0], a[0][0]*b[0][1]+a[0][1]*b[1][1]],
                     [a[1][0]*b[0][0]+a[1][1]*b[1][0], a[1][0]*b[0][1]+a[1][1]*b[1][1]]]
    
    def matrix_pow(a, n):
        if n == 1:
            return a
        elif (n % 2):
            return matrix_mul(a, matrix_pow(matrix_mul(a, a), n/2))
        else:
            return matrix_pow(matrix_mul(a, a), n/2)
    
    def fib_log(n):
        a = [[1, 1], [1,0]]
        return matrix_pow(a, n)[1][0]
    
  4. Graham said

    …and also apologies about overindenting the second line of matrix_mul. It didn’t look right in the comment box, but I should have left it alone.

    I really appreciate your excellent blog. Please keep up the good work!

  5. F. Carr said

    Here is a logarithmic version — but using base 3 instead of base 2.

  6. F. Carr said

    Oops, left off the link to codepad.org.

  7. Remco Niemeijer said

    Alternative (and, admittedly, slightly hacky, as evidenced by the warnings you get if you don’t specify -fno-warn-missing-methods and -fno-warn-orphans) Haskell solution for the logarithmic version, based on the fact that the default power operator is already logarithmic:

    {-# LANGUAGE FlexibleInstances #-}
    
    import Data.List
    
    instance Num a => Num [[a]] where
      a * b = [map (sum . zipWith (*) r) $ transpose b | r <- a]
    
    fiblog :: Int -> Integer
    fiblog n = ([[1,1],[1,0]] ^ n) !! 1 !! 0
    
  8. Christopher Oliver said

    An alternate scheme implementation with minor bit diddling.

    ;;; Computation of Fibonacci numbers with F(0) = 0 and F(1) = 1
    
    (define (naive-fib n)
      (if (< n 2) n (+ (naive-fib (- n 1)) (naive-fib (- n 2)))))
    
    (define (linear-fib n)
      (let loop ((a 0) (b 1) (n n))
        (if (< n 1) a (loop b (+ a b) (- n 1)))))
    
    ;; The following takes advantage of matrix symmetry and the constant
    ;; successor term to use cheaper steps for the square and multiply
    ;; operations.  Note that the matrix square takes Fib n to Fib 2n+1,
    ;; so the sense of even/odd in the power function is reversed from
    ;; the usual.
    (define (fast-fib n)
      (define (successor a b c)
        (list (+ a b)
    	  (+ b c)
    	  b))
    
      (define (square a b c)
        (let ((b-squared (* b b)))
          (list (+ (* a a) b-squared)
    	    (* b (+ a c))
    	    (+ b-squared (* c c)))))
    
      (define (matrix-fib n)
        (cond ((= n 1) '(2 1 1))
    	  ((even? n) (apply successor (matrix-fib (- n 1))))
    	  (else (apply square (matrix-fib (/ (- n 1) 2))))))
    
      (if (< n 1) 0 (caddr (matrix-fib n))))
    
  9. Christopher Oliver said

    And since the m[0,0] = m[1,1]+m[1,0]…

    (define (fast-fib n)
      (define (successor a b)
        (list (+ a b) a))
    
      (define (square a b) \
        (list (* a (+ a b b)) (+ (* a a) (* b b))))
    
      (define (matrix-fib n)
        (cond ((= n 0) '(1 0))
    	  ((even? n) (apply successor (matrix-fib (- n 1))))
    	  (else (apply square (matrix-fib (/ (- n 1) 2))))))
    
      (if (< n 1) 0 (car (matrix-fib (- n 1)))))
    

    That will teach me to proofread.

  10. Christopher Oliver said

    Err!!! delete that backslash. Praxis needs a way for responders to edit their mistakes!

  11. rahul said

    please give also solution in c langauge pleaseeeeeeeee.

Leave a Reply