Fibonacci Numbers
July 30, 2010
Here is the classic exponential algorithm, using recursion:
(define (fib-r n)
(cond ((zero? n) 0)
((< n 2) 1)
(else (+ (fib-r (- n 1)) (fib-r (- n 2))))))
And here is the classic linear algorithm, using iteration:
(define (fib-i n)
(let loop ((n n) (f-1 1) (f-2 0))
(if (zero? n) f-2
(loop (- n 1) (+ f-1 f-2) f-1))))
The logarithmic algorithm requires a powering function on matrices. We use the matrix datatype from the Standard Prelude and matrix multiplication from an earlier exercise to write this function that multiplies the matrix m by itself n times:
(define (matrix-power m n)
(cond ((= n 1) m)
((even? n)
(matrix-power
(matrix-multiply m m)
(/ n 2)))
(else (matrix-multiply m
(matrix-power
(matrix-multiply m m)
(/ (- n 1) 2))))))
Then the logarithmic-time fibonacci calculation is simple:
(define (fib-m n)
(if (zero? n) 0
(matrix-ref (matrix-power #(#(1 1) #(1 0)) n) 1 0)))
Here are some timings:
> (time (display (fib-r 25)) (newline))
75025
cpu time: 10 real time: 75 gc time: 0
> (time (display (fib-i 25000)) (newline))
2195438355...
cpu time: 80 real time: 518 gc time: 50
> (time (display (fib-m 25000)) (newline))
2195438355...
cpu time: 0 real time: 6 gc time: 0
Wow! You can run the program at http://programmingpraxis.codepad.org/YpIgShGQ.
Programming Praxis 
[…] Praxis – Fibonacci Numbers By Remco Niemeijer In today’s Programming Praxis we have to provide three different methods of calculating the ever-popular […]
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 !! 0I 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]…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!
Here is a logarithmic version — but using base 3 instead of base 2.
Oops, left off the link to codepad.org.
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 !! 0An 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))))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.
Err!!! delete that backslash. Praxis needs a way for responders to edit their mistakes!
please give also solution in c langauge pleaseeeeeeeee.
After writing this exercise, I came across Dijkstra’s paper describing a set of recurrence equations that calculate the nth fibonacci number. His method is orthogonal to the matrix method described above, but more convenient for computation. Note that Dijkstra starts the sequence differently; where we say F0=0, Dijkstra says F0=1. Here’s the code:
(define (fib n)(define (square x) (* x x))
(cond ((zero? n) 0) ((or (= n 1) (= n 2)) 1)
((even? n) (let* ((n2 (quotient n 2)) (n2-1 (- n2 1)))
(* (fib n2) (+ (* 2 (fib n2-1)) (fib n2)))))
(else (let* ((n2-1 (quotient n 2)) (n2 (+ n2-1 1)))
(+ (square (fib n2-1)) (square (fib n2)))))))
I’ve seen difference equations (e.g., fibonacci) expressed as matrices before, but this was the first time I’ve seen the concepts of logarithmic-time exponentiation, matrices, and difference equations all in one problem. This was a genuinely entertaining exercise; thanks for posting.
Since the core of this problem was essentially finding a O(log(n)) algorithm for matrix multiplication, I also wrote a O(1) solution based upon diagonalization. Looks like it works, but for very large values, it may suffer from accuracy issues due to the use of inexact numbers.
(define (transpose A) (apply map list A)) ;; crude matrix-mult for 2x2 (define (matrix-mult A B) (let ((trans-B (transpose B))) (cons (list (apply + (map * (car A) (car trans-B))) (apply + (map * (car A) (cadr trans-B)))) (list (list (apply + (map * (cadr A) (car trans-B))) (apply + (map * (cadr A) (cadr trans-B))))) ))) (define (log_fib n) (cond ((>= 1 n) '((1 1) (1 0))) ((even? n) (let ((fib-matrix (log_fib (/ n 2)))) (matrix-mult fib-matrix fib-matrix))) (else (let ((fib-matrix (log_fib (quotient n 2)))) (matrix-mult (matrix-mult fib-matrix fib-matrix) '((1 1) (1 0))))))) (define (PDP_fib n) (let* ((eig1 (* 1/2 (+ 1 (sqrt 5)))) (eig2 (* 1/2 (- 1 (sqrt 5)))) (P (list (list eig1 eig2) (list 1 1))) (D (list (list (expt eig1 n) 0) (list 0 (expt eig2 n )))) (P^-1 (let ((det (/ 1 (- eig1 eig2)))) (list (list det (* det eig2 -1)) (list (* det -1) (* det eig1)))))) (round (cadr (car (matrix-mult (matrix-mult P D) P^-1))))))[ I fixed the formatting. In the future, look at the code formatting how-to for instructions. Thanks for your kind words, and for an interesting variation on the problem. — PP ]
public void fibonacci(int firstTerm, int secTerm)
{
if (nextTerm == 75025) return;
nextTerm = firstTerm + secTerm;
System.out.println(counter+”. “+nextTerm);
fibonacci(secTerm, nextTerm);
}
Goes to 25th term, based on firstTerm = 0 and secTerm = 1.
Gil Martinez, I think your O(1) approximation is about as honest as the “O(log n)” algorithm we see here. It’s actually impossible to achieve O(log n) time for calculating arbitrarily large Fibonacci numbers: to represent fib(n) requires O(log(fib(n))) bits. But as shown at this site, fib(n)=round(φn / √5), which is Θ(φn), so we need Θ(log(φn))=Θ(n) bits, which we can’t possibly produce in O(log n) time. What I don’t know is whether the naive algorithm is actually O(n) or O(n log n).
Those were all supposed to be φ^n. Not sure where that went screwy.
I missed something more important: The PDPFib function is defined using (expt eig1 n). I don’t see how that function can possibly be O(1), even assuming that addition, subtraction, multiplication and division are.