Compute Fibonacci numbers FAST!

The Fibonacci numbers are a commonly-given example in Computer Science and in Mathematics. They will undoubtably make an appearance in every single introductory programming course. It is therefore perhaps surprising that an efficient method for computing them is so little-known! Instead, we make do with methods which take time O^~(n) (or even, horror of horrors, time O^~(F_n)!) to compute F_n.

The recursive definition

F₀=0
F₁=1
F_n+1=F_n+F_n-1

leads to a straightforward implementation. Unfortunately, this implementation is horribly inefficient: to compute F_n+1, you need to compute F_n and F_n-1, and perform an extra addition. Solving the recurrence relation yields time O^~(F_n) to compute F_n -- this is time exponential in n, therefore exponentially exponential in the size of the input n.

Most courses will therefore teach a somewhat faster algorithm: at any given step, keep variables x and y, with x being F_n and y being F_n+1. The loop updates by saying

x' = y
y' = x+y

(a temporary variable will often be needed to prevent one update clobbering the other). This loop requires about n iterations to compute F_n, so the time required is O^~(n) -- linear in n, but exponential in the size of n.

But this algorithm too is rubbish! Binet's formula gives us a hint: by computing powers of some irrational quadratic surd, we can get the answer. And lots of people have seen how to compute the n'th power in time O^~(log n). Why not apply it here, and be done in logarithmic time?

Unfortunately, this approach cannot work: you'll need to compute an indecently large number of digits of sqrt(5) (or alternatively of the golden ratio φ), which will nullify your advantages. We need to take powers of something else.

That something is the square 2x2 matrix

[[1 1]
 [1 0]]

This matrix (call it G) has the property that Gⁿ is exactly

[[Fn+1 F n ]
 [F n  Fn-1]],

as can easily proved by induction. (People with a background in linear algebra might note that the matrix's characteristic polynomial is the familiar x²-x-1, so that its eigenvalues are precisely the golden ratio φ and its complement -φ+1/2. Really bored people with a background in linear algebra will go on to diagonalize the matrix and derive Binet's formula. Really really bored people will go on to derive the general solution for a linear recurrence from this technique.)

So, to compute F_n in time O^~(log n), we need to show how to compute Gⁿ with O(log n) matrix multiplications. But that's easy (see also fast exponentiation for an explanation of this process)!

G⁰ = I
G^2k+1 = G^2k*G
G^2k = (G^k)²

So we can compute Gⁿ efficiently; the rest is merely coding.

This is probably about the most efficient method of calculating Fibonacci numbers. It's exponentially better than "efficient" linear algorithm, and exponentially exponentially better than the naïve translation of the recurrence relation!

---

Note on "O^~(f(n))"

This formalism, used above, will not be too familiar. It means that O(f(n)) arithmetic operations on words are required, where words are assumed to be "big enough" to hold the result. Of course this isn't true -- words have some limited size. If we wish to have unbounded n, we must replace these operations with operations on variable-width numbers: a 256-bit adder will be slower than a 32-bit adder (even if a carry lookahead adder is used); in a program, no parallelism is possible, and the slowdown is even greater. This is particularly important with number which grow exponentially with n, like the Fibonacci numbers: the number of digits required for F_n is O(n).

Still, counting arithmetic operations is easy to do, and is accurate enough for most purposes. If required, we can translate the time bounds above to absolute bounds ("time O(g(n))") by multiplying by appropriate factors. The logarithmic algorithm remains much better.

If you find the matrix implementation above intimidating, here's how to get the same O(log n) behavior* using an amazing relationship between Fibonacci numbers and the Golden Ratio. The Golden Ratio, typically represented by the Greek letter phi (Φ), is defined as

Φ = (1 + 5^1/2) / 2,

or approximately 1.618. We'll also use the conjugate of Φ, which we'll call Φ^:

Φ^ = (1 - 5^1/2) / 2,

or approximately -.618. The formula for the nth Fibonacci number, then, is simply:

F_n = (Φ^ⁿ - Φⁿ) / 5^1/2

We can further simplify the formula by getting rid of the conjugate, since its contribution is always pretty negligible. Because |Φ^| < 1, we know |Φ^ⁿ| / 5^1/2 < 1/5^1/2 < 1/2, so the subtracted portion will always be less than a half. Therefore we can say F_n = Φⁿ / 5^1/2, rounded to the nearest integer. In other words, to get the nth Fibonacci number, all you have to do is take the Golden Ratio to the nth power, divide it by the square root of 5, and round it off. Amazing, right?

But why is any of this the case? Well, Fibonacci numbers and the Golden Ratio both occur with freakish regularity in nature, so it makes sense that they would be cosmically linked somehow. Ok, so the real answer is I have no freaking clue why**. But logic shows us things we can't necessarily grasp intuitively, and I can at least prove that it's true.

* The original version of this writeup made the somewhat smug assertion that we could get O(1) time from this "algorithm", but pfft set me straight: computing the nth power of phi still requires log n time. This is why I shouldn't do technical nodes.

** For a real explanation of the relationship between the Fibonacci numbers and the Golden Ratio, see mblase's definitive writeup under the latter.

---

Proof

The recursive definition of the Fibonacci sequence is:

F₀ = 0
F₁ = 1
F_n+2 = F_n + F_n+1 (n >= 0)

Let P(n) be the proposition that F_n = (Φⁿ - Φ^ⁿ) / 5^1/2.

For n of 1 and 2: F₁ is defined to be 1. (Φ¹ - Φ^¹) / 5^1/2 = 1, so P(1) is true. F₂ is defined to be F₀ + F₁ = 0 + 1 = 1. (Φ² - Φ^²) / 5^1/2 = 1, so P(2) is true as well.

Now for the inductive step. Assume P(k) and P(k+1) to be true; that is, F_k = (Φ^k - Φ^^k) / 5^1/2 and F_k+1 = (Φ^k+1 - Φ^^k+1) / 5^1/2. We will show that P(k+2) can be inferred from this. By definition, F_k+2 = F_k + F_k+1. Based on the above assumption, we can substitute this for the right side of the equation:

F_k+2 = (Φ^k - Φ^^k) / 5^1/2 + (Φ^k+1 - Φ^^k+1) / 5^1/2

With some algebraic manipulation,

F_k+2 = ((Φ^k - Φ^^k) + (Φ^k+1 - Φ^^k+1)) / 5^1/2
F_k+2 = (Φ^k + Φ^k+1 - Φ^^k - Φ^^k+1) / 5^1/2
F_k+2 = (Φ^k(1 + Φ) - Φ^^k(1 + Φ^) / 5^1/2
F_k+2 = (Φ^kΦ² - Φ^^kΦ^² / 5^1/2 ***
F_k+2 = (Φ^k+2 - Φ^^k+2) / 5^1/2

Which is exactly what we were proposing to show. Therefore (P(k) & P(k+1)) --> P(k+2) and the proof is complete.

*** Show that Φ² = 1 + Φ:

Φ² = ((1 + 5^1/2)/2)²
Φ² = (1 + 2*5^1/2 + 5) / 4
Φ² = (6 + 2*5^1/2) / 4
Φ² = (3 + 5^1/2) / 2
Φ² = (2 + 1 + 5^1/2) / 2
Φ² = 2/2 + (1 + 5^1/2) / 2
Φ² = 1 + Φ

(The same process can be used to show that Φ^² = 1 + Φ^)

---

ariels says "One problem with using Binet's formula (as you do) is that in order to perform accurate calculations of φ^n/sqrt(5) you need to expand the two quadratic surds (φ and sqrt(5)) to increasing precisions. I think you run out of precision quite quickly with floats or even doubles..."

I've tried all the methods above, but the method used by the GNU Multiple Precision (GMP) arithmetic library is at least five times faster than its nearest competitor in Ruby. Because GMP is optimized for speed and written in C, I assume the method works quickly in other languages as well.

To find the n^th Fibonacci number, use the following recursive definition :

• F₀ = 0
• F₁ = 1
• When n is an even number :
     F_n = F_{n / 2}(F_{n / 2} + 2 F_{(n−2) / 2})
• When n, modulo 4, is 1 :
     F_n = (2 F_{(n−1) / 2} + F_{(n−3) / 2})(2 F_{(n−1) / 2} − F_{(n−3) / 2}) + 2
• When n, modulo 4, is 3 :
     F_n = (2 F_{(n−1) / 2} + F_{(n−3) / 2})(2 F_{(n−1) / 2} − F_{(n−3) / 2}) − 2

---

Here's the Ruby code (using a hash table for memoization) :

class Integer
  FibonacciComputer = Hash.new do |fibonacci, n|
    if fibonacci.has_key?(n - 1) and fibonacci.has_key?(n - 2)
      fibonacci[n] = fibonacci[n - 1] + fibonacci[n - 2]
    elsif fibonacci.has_key?(n + 1) and fibonacci.has_key?(n + 2)
      fibonacci[n] = fibonacci[n + 2] - fibonacci[n + 1]
    else
      half_n = n.div(2)
      case n.modulo(4)
        when 1
          fibonacci[n] = (2*fibonacci[half_n] + fibonacci[half_n - 1])*(2*fibonacci[half_n] - fibonacci[half_n - 1]) + 2
        when 3
          fibonacci[n] = (2*fibonacci[half_n] + fibonacci[half_n - 1])*(2*fibonacci[half_n] - fibonacci[half_n - 1]) - 2
        else
          fibonacci[n] = fibonacci[half_n]*(fibonacci[half_n] + 2*fibonacci[half_n - 1])
      end
    end
  end
  FibonacciComputer[0] = 0
  FibonacciComputer[1] = 1
  FibonacciComputer.freeze
  
  def fib
    return FibonacciComputer.dup[self]
  end
end

puts "The 1,000th Fibonacci number is #{ 1_000.fib }."

Try it out at http://www.ruby.ch/interpreter/rubyinterpreter.shtml.