Fibonacci Numbers in the Real World

Introduction

I learned how to program computers from a utilitarian point of view: fortran¹ for fast numerics, Python for exploring data and making website back ends, some javascript for the front end.

I’d heard about the famed power and beauty of the lisp family of languages, and even learned how they worked a bit, out of curiosity. (A similar curiosity drive me to learn the basics of Haskell). I’d read the inspirational essays, and been intrigued, but two things finally led me to begin learning a lisp language in earnest: the growth of Clojure, a uniquely modern, “practical” lisp, and this:

[1]

(def fibs 
   (lazy-cat [0 1] (map + fibs (rest fibs))))

This is a line of Clojure that generates the Fibonacci sequence. If you execute this definition in your REPL (interactive Clojure prompt) and then type fibs, you will see the infinite Fibonacci sequence begin to print out. (I usually do (set! print-length 200) in the REPL to limit how much of a big list gets dumped on the screen.) To see just the single Fibonacci number that you are interested in, you might want to say (nth fibs 100)

This is an example of corecursion. I had seen and used recursion before, in its usual form of a function or subroutine that invokes itself. But here was a new and strange species of recursion that did not explicitly use any type of callable at all; that simply defined a list and caused computation to happen just by saying its name. The idea that a calculation could be described as the recursive generation of a sequence without the use of functions was, to me, a new and fascinating idea.

But beauty, as we know, is not always practical. As pointed out here, line [1] consumes too much memory. A more practical, if less beautiful, translation of this idea is to define an equivalent function:

[2]

(defn fibn [] 
   (map first (iterate 
              (fn [[a b]] [b (+ a b)]) [0 1])))

Now the \(n^{\rm th}\) Fibonacci number is extracted with (last (take (inc n) (fibn)))).

A Wakeup Call?

The reverie of the next few days, exploring Clojure and other lisp hackers’ elegant solutions to the generation of Fibonacci numbers and other problems, was interrupted by a bracing article by Evan Miller entitled “The Mathematical Hacker.” I highly recommend this elegantly written essay, especially if you’ve found what I’ve written here so far to be at all interesting. I recommend it even more urgently if you’ve found the above not very interesting because you’ve seen it all already.

One of the chief points of Miller’s article is that lisp hackers are perhaps a bit too fond of their ingenious solutions to problems such as the generation of Fibonacci numbers. After all, the clever recursive algorithms typically run in time roughly linear in \(n\), and only make sense for integer \(n\).

As Miller points out, the problem of finding any Fibonacci number you want yielded to an elegant, closed-form solution in the 18^th (he has 17^th) century. The formula, and related results, are derived in many places; I like the treatment in Knuth, which is characteristically lucid and readable. There we have, on page 82,

[3] \[F_n = \frac{1}{\sqrt{5}} (φ^n - \hat{φ}^n),\]

where \(F_n\) is the \(n^{\rm th}\) Fibonacci number, \(φ\) is the golden ratio \(\frac{1}{2}(1 + \sqrt{5})\) and \(\hat{φ} = \frac{1}{2}(1 - \sqrt{5})\). This is just as beautiful as [1], and, in its way, just as surprising; the treatment in Knuth’s book makes it clear where the mysterious relationship between this numerical sequence and the \(\sqrt{5}\) comes from.

Miller presents a C program, amounting to a direct translation of [3], that he claims is an obviously better way to produce Fibonacci numbers, for several reasons. The direct use of [3] would seem to run in constant time, regardless of the value of \(n\), because it executes a fixed number of arithmetic operations; there is no looping or recursion that increases the number of operations as \(n\) increases. In contrast, clever lisp-hacker type solutions such as [1] typically take time linear in \(n\): it will take roughly twice as long to get the 100^th Fibonacci number as the 50^th.

Another reason using the analytical formula is preferable, according to Miller, is that it encourages us to think about the mathematics underlying it, and wonder where the \(\sqrt{5}\) comes from, or what it means to plug in non-integer values for \(n\). We get involved in the rich history and culture of mathematics, and drawn out of the “programmer’s cubicle,” a supposedly intellectually impoverished place in which Miller claims the lisp hacker spends too much time.

Not So Fast

Miller may be right when he says that we need an “infusion of applied mathematics into hacker education.” It is probably true that many programmers would be able to reach better solutions for some of their problems if they stopped thinking of them as programming problems, and were able to draw upon a fuller understanding of our mathematical heritage. My personal prejudice is that everyone would be better off for this, from painters to biologists, so I am an easy sell for his underlying message. However, Miller seems to have chosen an example that undercuts his agenda somewhat. Further examination shows that the lisp hackers may have the last laugh in this case, and most likely in a large class of related cases.

First, let’s look at the mathematical content itself. Equation [3] is surprising, beautiful, and suggestive, it is true. But the mathematics required to derive it is fairly elementary; if you are comfortable with manipulating infinite series, you can follow it. In other words, anyone who has completed a beginning calculus (or even what is called in some places “pre-calculus”) class will encounter no fundamentally new mathematics in the derivation of the formula, although the particular techniques of proof might seem novel, if the student has not seen generating functions before.

The line of code in [1] is also, at least in my opinion, surprising, beautiful, and suggestive. Further, my interest in understanding how it works, and later, the ideas behind it, led me into unfamiliar realms of computer science, which is a branch of mathematics itself. I say unfamiliar, but that is relative to my background in physics and mathematics; the ideas behind corecursion and related concepts may be old hat to those, unlike me, with a solid and fairly recent grounding in computer science. So I’d like to turn Miller’s suggestion around and suggest that learning how lisp hackers approach these problems has led me out of my physicist’s cubicle with all its applied math and fortran programs and down a road strewn with new and powerful tools. Perhaps, as well as an infusion of applied math into hacker education, we might benefit from a carefully metered IV drip of hacker knowledge into science and engineering education.

Aside from these considerations, Miller’s example might not be the best choice for illustrating his points. If we run the C program that he recommends for calculating Fibonacci numbers, we find that it works great, using constant time, and not much of it, to generate the sequence up to \(F_{46} = 1836311903\). After this, it falls apart, because this is as far as you can go with the long int datatype used in the program.² You might be able to get larger numbers using the long long int datatype in some environments, but to go arbitrarily far out in the sequence you would need to resort to an arbitrary precision numerical library.

Since arithmetic operations on arbitrary precision numbers take time that increases with \(N\), the number of digits in the operands, a method that uses a formula such as [3] can not take time constant in \(N\), unless we limit ourselves, in this case, to \(F_{46}\). But if we were merely interested in small Fibonacci numbers, the complexity of our algorithm would hardly matter, unless it were really extremely inefficient. Therefore one of the advantages of the method recommended by Miller evaporates, or at least is not as obvious as he claims.

There is one sense in which a direct application of [3] can be said to produce Fibonacci numbers in constant time, if we relax what we mean by a Fibonacci number. If we use, instead of integers, floating point numbers, then our results will take the form of a number with a certain fixed precision and a magnitude given by the exponent. How much precision, and the range of exponents, depends on whether we decide to use single precision or double precision floats. So if we only care about the first handful of significant digits and the general magnitude of our Fibonacci numbers, rather than their exact representation, Miller’s recommended method will produce these answers in constant time. But this is not what we usually think of when we talk about generating Fibonacci numbers, factorials, etc.

Assuming from here on that our problem is to calculate exact, large Fibonacci numbers, it is not immediately clear whether a direct use of [3], or an iterative algorithm like [2],³ is faster for large \(n\). Although the “mathematical” solution is pure arithmetic and has no looping, it does involve exponentiation, which may become rapidly more expensive as \(n\) (and therefore \(N\)) increases. Solution [2], although it requires more operations as \(n\) increases, only involves addition, which is at worst \(\mathcal O(N)\), inserting numbers into a vector, and function application. In other words, a naive complexity analysis of the mathematics based solution might conclude that it executes in constant time, but some of the complexity is hidden within the arithmetic operations themselves, when they operate on arbitrary precision operands.

Experiments

Now let’s crank out some Fibonacci numbers. We would like to find out how the applied mathematician’s approach compares with the lisp hacker’s approach. In order for the comparison of execution times to be meaningful, we’ll code both solutions in Clojure. This might be a good place to reiterate that I am a rank Clojure beginner. These bits of code will do what they are supposed to do, but might be stupid in any number of ways about which I have no clue.

If you are following along you will need to include these dependencies in your project file:

[org.apfloat/apfloat "1.6.3"]
[org.clojure/math.numeric-tower "0.0.2"]

And, in your REPL or source file, put

(require '[clojure.math.numeric-tower :as math])
(import '(org.apfloat Apfloat ApfloatMath))

The numeric-tower library gives us exponentiation and rounding functions. The apfloat import gives us access to Java arbitrary precision numeric classes. We’ll find out below why we need to use these at all, even though Clojure has its own arbitrary precision numbers.

For timing computations I used the well regarded Criterium benchmarking library. This produces more consistent and meaningful results than simply using the time function. If you want to use Criterium, then add [criterium “0.4.1”] to the list of dependencies in your project file, and include (use ’criterium.core) or the equivalent in your source file. For all the timings reported here I used Criterium’s quick-bench function and ran the code on an Intel 1.83 GHz dual-core processor with one CPU core turned off.

Experimenting with the definition of fibn in [2] above, we discover that the Fibonacci numbers are accurately calculated up to \(F_{91}\), above which Clojure has the good sense to say ArithmeticException integer overflow rather than producing meaningless results. This is about what you would expect for 64 bit integers.

If we replace integers in the [0 1] vector initialization in [2] with floats ([0.0 1.0]), we can get a numerical output up to \(F_{1476}\) = 1.3069892237633987 × 10³⁰⁸, above which the output is “Infinity.” This is close to the largest representable 64 bit float. But as pointed out above, this is not what we want. Indeed, if we only need a rough idea of the magnitude of the Fibonacci number, there are simple enough analytic ways to estimate that.

Figure 1 gives the time to calculate a Fibonacci number vs the number in question, \(F_n\), using the float version of fibn. The inset graph shows the integer version of the same algorithm. By coincidence the time to calculate a Fibonacci number using either version, expressed in microseconds, is almost exactly the same as the Fibonacci number itself. The iterative approach is clearly \(\mathcal O(N)\), as stated by Miller in his article. The timing results are almost perfectly linear, as can be seen by the straight line fit plotted on top of the results.

As an aside, a direct Clojure translation of [3] using native datatypes returns accurate Fibonacci numbers up to \(F_{70}\) in a constant time of 1.7 μs.

In order to extend our calculation of exact Fibonacci numbers past \(F_{91}\), we need to use datatypes with greater precision than that available from long ints, as alluded to in the previous section. We’ll make a small change to the definition of fibn in [2]:

[4]

(defn fibo [] 
   (map first (iterate 
              (fn [[a b]] [b (+ a b)]) [0N 1N])))

We’ve used the 0N and 01 notations for the initialization of the array; this is Clojure’s notation for BigInts, arbitrary precision integers. Now, as the array is extended, subsequent elements will be promoted to BigInts, and our Fibonacci numbers will have all the digits that they need. As before, if we say

[5]

(defn ffibo [n]
   (last (take (inc n) (fibo))))

then (ffibo n) will return the \(n\)^th Fibonacci number.

This function now returns any Fibonacci number we want, perfectly accurately, happily spewing thousands of digits onto the screen. Figure 2 shows fibo’s execution timings. Now the run time as a function of Fibonacci number has a small quadratic component, where the ratio of b to c in the fitting curve is 12,052. The origin of this quadratic growth is the convolution of the linear dependence arising from the iteration with the, presumably, \(\mathcal O(N)\) dependence of addition of the arbitrary precision numbers containing \(N\) digits.

So much for the lisp hacker solution. A straightforward translation of [3] into Clojure would be

(defn afib [n]
   (let [sq5 (Math/pow 5 0.5M), rsq5 (/ 1 sq5)]
   (math/round (* rsq5 
   (- (Math/pow (/ (+ 1 sq5) 2) n) 
      (Math/pow (/ (- 1 sq5) 2) n)))))) 

The notation 0.5M means 0.5 is a BigDecimal, which is another unlimited precision Clojure datatype.

Experimenting with afib shows that it produces the correct result up to (afib 70), but for Fibonacci numbers higher than this, the least significant digits start to go wrong. Further research shows that precision is lost by the Math/pow function, which returns a double even when operating on arbitrary precision operands. The math/expt function, another option, will return a BigDecimal when raising a BigDecimal to an integer power, but in all other cases returns a double. Therefore we could retain arbitrary precision in afib if we could precompute \(\sqrt{5}\) to the required precision for each Fibonacci number to be calculated and store it in a BigDecimal. I was pointed to another way to handle these problems by Takahiro Noda, in answer to the only question I’ve ever asked on Stack Overflow. There is an arbitrary precision arithmetic library in the Java world that, of course, we can use from Clojure. This is what the apfloat dependencies, that we listed at the beginning of this section, are for. Translating the definition of afib into a form that uses the apfloat library, I come up with

(defn tptfib [n]
   (let [pcn (+ 10 (int (* 0.21 n)))
         one (Apfloat. 1M pcn) 
         two (Apfloat. 2M pcn) 
         asq5 (ApfloatMath/pow (Apfloat. 5M pcn) (Apfloat. 0.5M pcn))] 
   (defn tfib [n]
      (.ceil (.multiply (.divide one asq5)
      (.subtract (ApfloatMath/pow (.divide (.add one asq5) two) n) 
      (ApfloatMath/pow (.divide (.subtract one asq5) two) n))))))) 

This is a function that, when applied to an integer \(n\), returns another function, tfib, that returns the \(n\)^th Fibonacci number with the required number of digits. From [3] we can see that, once \(n\) gets moderately large, \(F_n\) has about \(0.21 n\) digits (when written as a number in base 10); in the definition of tptfib we’ve added 10 additional digits to take care of the extra precision required by the intermediate expressions.

The code in [2] works, producing giant Fibonacci numbers just as readily as our fibo function. Figure 3 compares its timing with fibo, showing a more chaotic dependence on \(n\) but generally following the same quadratic trend. The “lisp hacker” solution performs better throughout most of the range up to about \(F_n\) = 25,000, often producing the result in half the computation time.

Note that in Figure 3 we are actually plotting the execution times of the function tfib returned by tptfib, so that we are not including the time required to instantiate the Apfloat constants. It was with timing studies in mind that tptfib was written this way.

The chaotic nature of tptfib’s execution times, which is largely repeatable across timing experiments using both the Criterium library and the time function, may be a symptom of the complicated algorithms used in the apfloat library. By “repeatable” I mean that the curves do not vary much if the timing experiment is repeated, revealing a genuinely complicated dependence of apfloat arithmetic execution times on the sizes of the operands.

The most obvious remaining possible experiment is to implement the algorithm in [2] using apfloat numbers rather than Clojure’s BigInts or BigDecimals. But I can see no practical way to do this, as new numbers with increasing precision would need to be continually instantiated as the list was built up. Alternatively, the list could be initialized with sufficient precision to hold the final result of interest, using, necessarily, more precision than required for each previous member of the list. This I have tried, and discovered a similar quadratic dependence as in the timings of fibo, but that the time required was roughly an order of magnitude greater. Therefore this approach is a non-starter. Finally, it should be mentioned that if your goal is to calculate Fibonacci numbers as efficiently as possible there are other well-known algorithms that are probably better choices than any looked at here.

Summary

I agree enthusiastically with all that Evan Miller says in his article about the importance of mathematics and the desirability of mathematical education for hackers. When you face the details of particular computations, however, you need to grapple with the gritty facts that the numbers inside a computer are not the numbers inside our heads (or in the Platonic realm), symbolized in the formulas of the annals of mathematics. Only a small, finite set of numbers have an exact representation as native datatypes inside the computer, for which there are efficient machine instructions to perform arithmetic upon them. Arithmetic on large integers, as we’ve seen above, takes time that grows with the size of the operands, as they must be broken into smaller parts or taken digit by digit. Any analysis that concludes that an algorithm takes constant time in \(n\) to calculate the \(n^{\rm th}\) Fibonacci number is almost certainly incorrect, therefore, as it considers a multiplication, for example, to take a fixed time for any two operands. And for this reason, sometimes the clever algorithms devised by the lisp hackers actually work better than a direct application of the elegant mathematical formulas that would appear, at first glance, to be better solutions. These clever “lisp hacker” solutions, furthermore, can be beautiful and provocative exemplars in their own right of another branch of the great, eternal, and eternally developing tree of mathematics.

Further Reading

Note added April 26^th, 2016: I just discovered the surprising and delightful article on “Paul’s Blog” called “An integer formula for Fibonacci numbers.” If you found this at all interesting, you might want to have a look at that.

See also (added July 26^th, 2019): Finding Fibonacci numbers using Linear Algebra

Added April 23^rd, 2020: This is quite germane and very nifty.

Appeared Jul. 30^th, 2023: Good explanations of the matrix methods for calculating Fibonacci numbers, with interesting animations.

From September, 2024: Sequences whose ratios converge to metallic constants