Real World Benchmarking
by Jay Coskey, 2011-09-12

The Fibonacci sequence has a simple definition: F(0) = 0; F(1) = 1; F(n) = F(n – 1) + F(n -2).  A naive translation from the definition to an algorithm for computing individual elements performs poorly, due to heavy branching.  Without memoization, this algorithm has time complexity Omega(phi^n), where phi = (1+sqrt(5))/2.  (See [1].)  With memoization, that can be dropped to linear time.

With such a simple definition, isn't there a constant-time formula for F(n)?  Yes, strictly speaking, there is.  By viewing the problem as iterated 2x2 matrix multiplication, one can use the eigenvalues and eigenvectors of the associated matrix to determine a closed formula for F(n).  This formula is called Binet's formula:
	F(n) = [phi^n - (-phi)^(-n)] / sqrt(5).
(A simplification involves just computing round(phi^n / sqrt(5)).  Using Binet's formula requires the use of either floating-point arithmetic, which causes errors for large values of n due to rounding, or so-called infinite precision arithmetic, which is significantly slower than methods that employ integer arithmetic.

By using integers and finding ways to reduce evaluation steps, we can achieve O(log n)  performance.  Here we compare the performance of three algorithms:
	1.	Extension.  Perform arithmetic in an extension (i.e., a superset) of the integers that contains the square root of five.  This extension of the integers is written as Z[sqrt(5)], and called the extension ring of the integers generated by sqrt(5) (See [2].)  With this extension ring, we can manipulate square roots of 5 while still dealing only with the integer arithmetic.
			* See the file Fibonacci/Extended.hs.
	2.	Matrix.  Use the iterated matrix multiplication mentioned earlier.  This reduces the computation of F(n) to a simple matrix exponent problem.
			* See the file Fibonacci/Matrix.hs.
	3.	Recurrence.  Use a recurrence relation that expresses F(n) in terms of series with much lower indices.  (See [3].)
			* See the file Fibonacci/Recurrence.hs.

The programming language Haskell lends itself well to this problem, given its built-in support for multi-precision integer arithmetic, using the Integer type.  Since Haskell uses lazy evaluation, care must be taken in Haskell to ensure that computation time span measurements capture the time of the entire computation.  The Criterion benchmarking library for Haskell accomplishes this by allowing the caller of a benchmark to specify whether the result should be in normal form (NF) or weak head normal form (WHNF).  (See [4].)

As seen in the performance histogram, FibonacciPerfHistogram.png, all three methods have roughly the same time complexity, though the method based on matrix multiplication (#2), is roughly 2x to 5x faster for larger values of n.

The attached data set has no data point for the Extension method for n = 1 billion.  This is because an "out of memory error" was raised before the result could be obtained [5].  The algorithm does not branch excessively, since the exponentiation operation (***) takes care of one or two bits of the exponent per call.  But the Integer values being handled are on the order of 2 ^x, where x = log_2(2 * phi) > 1.  Operations on numbers with close to one billion bits are resource-intensive.  But keep in mind that that's quite far beyond the performance limits of the heavily branching naive algorithm, which becomes intractable at around n = 50.

[1] MIT's Introduction to Algorithms, Lecture 3: Divide and Conquer.  Viewable as a video at http://www.catonmat.net/blog/mit-introduction-to-algorithms-part-two/

[2] http://mathworld.wolfram.com/ExtensionRing.html 

[3] http://en.wikipedia.org/wiki/Fibonacci_number#Identities 

[4] http://stackoverflow.com/questions/6872898/haskell-what-is-weak-head-normal-form/6889335#6889335 

[5] These computations were performed on a 32-bit computer with 2 Dual-Core 3.40GHz Pentium CPUs and 2 GB RAM.