Several people suggested that Binet’s closed-form formula for Fibonacci numbers might lead to an even faster algorithm. That’s an interesting idea, which we’re going to explore in this post. So, what is Binet’s formula? It goes like this:

where is the golden ratio.
I’d like to explain why this formula works, because it’s a lovely story; and, despite what you might think from most of the explanations on the web, it’s perfectly possible to understand and appreciate it without knowing what “eigenvalue” means. But I fear that would be too long a digression here, so I’ll save it for another post.

There is a nice trick we can use to simplify the computation of this expression, taking advantage of the fact we know it always comes out as an integer. Since φ is about 1.618, its negative cousin (1-φ) is about -0.618, which means that the term

is always between -1/2 and 1/2. So we can just compute

and round to the nearest integer.

The simplest way to implement this idea is to use ordinary floating-point arithmetic. For example, in C we might write:

This is simple and fast. The only trouble is that it almost always gives the wrong answer. On my machine the answer is wrong whenever n > 70. For example fib_fixed(71) gives 308061521170130, whereas the correct answer is 308061521170129, and it only gets worse from there.

We can improve the situation a little by tweaking the computation, but, however cunningly we do that, we’ll soon run into a fundamental obstacle: fib(79) is 14472334024676221, which can’t even be represented in an IEEE double. You can see that illustrated as a comedy Google calculator fail.

So if we want to implement this in a general way then we need to use arbitrary-precision floating-point arithmetic. For some reason this has not gone mainstream in the way that big ints have done, but the GMP library has a good implementation, so we’ll use that. Here’s the code:

This works fine, gives correct answers, and is reasonably fast. But how does it compare to the integer methods? For comparison, I rewrote the fib_fast algorithm from my last post to use GMP integers. Click to see the code:

Then I made a graph showing the execution time of both methods, against the input n. The vertical axis shows elapsed time, in ticks, when running on my computer. (But it doesn’t really matter what it is, since we are interested in the relationship between the speed of the two methods, rather than the exact speed of either.)

Clearly the integer method is faster, but the curves are not obviously different shapes. A log-log plot shows the relationship more clearly:

which supports the idea (consistent with what you’d expect in theory) that both algorithms have the same big-O complexity, but the integer method is consistently about five times faster.

Floating-point arithmetic has a reputation for being difficult and unreliable; a technique that is good for finding approximate answers quickly, but best avoided when a precise answer is required. And actually I think that’s a fair assessment on the whole, but it’s interesting to see how arbitrary-precision floating point arithmetic can give us an exact answer very easily in this case. Of course we lose a lot of the speed advantage that floating-point conventionally offers, because our hardware doesn’t have the same fast native support for these “big floats” that it does for IEEE doubles.

Another approach

A related idea is to take advantage of the fact that all the calculations produce numbers of a special form, e.g. a+bφ for integers a and b. This works essentially because φ2 = 1 + φ, and so any polynomial in φ can be reduced to that form.

This turns out to be precisely equivalent to the matrix technique we used last time, because φn = fib(n-1) + fib(n)φ, and so computing powers of φ in this representation is just directly equivalent to computing Fibonacci numbers.

A variation of this – as suggested on HN – is to use the fact that φn is always of the form

This article made me check my Fibonacci routines in Delphi, and sure enough, they were slightly off at “high” values like 71 and 79. However, Delphi and its ancestors have access to Intel’s extended precision reals (extended type; most likely “long double” for most C compilers). Both the integer and floating point functions were using extendeds (long doubles) where appropriate, but the constants I had defined for Phi and SqrtOf5 were only defined to the precision of a double. After looking up the current definitions of those values (in Wikipedia, which changed a lot since those routines were written) and adding a few more digits of precision, the old routines provide correct values over a greater range of their domains. It’s not important enough for me to test thoroughly right now, and I’m short on time… but it looks pretty good on quick tests.

Oh yeah… I don’t return doubles; either int64’s (long long) or extended (long double). The long long return type doesn’t overflow until the input paramter is greater than 92! Of course the extended function will never show more than 20 digits of precision… One of these days I’ll find a good arbitrary-precision lib for floats for old Delphis…

I think there’s a small problem in your code for binet’s formula. You estimate the value of the base 2 log of the golden ratio right, but when you multiply by n and it gets typecast to an unsigned long (per the gmp docs, mp_bitcnt_t is unsigned long) many (if not all) C compilers floor the result to typecast. You can loose some precision because you’re just shy of what you need. This was evidenced by some large (circa 10^6) computations compared against an iterative integer computation (unrolling the tail recursion of a recursive solution) also using gmp.