The Nth Fibonacci Number in O(log N)

[cperciva]

It's easy to compute things quickly if you pretend that you can do arbitary-precision arithmetic in constant time. Want a quadrillion digits of Pi? Easy, just use the Brent–Salamin AGM algorithm and you'll get them in about 500 arithmetic operations. NP-complete problems? No problem at all, rephrase them as subset-sum problems and use the polynomial-integer-operations dynamic programming algorithm. I don't think you can solve PSPACE in a polynomial number of integer operations, but I can't see how to prove it immediately.

In "real" terms, computing the Nth fibonacci number takes O(M(N)) = O(N log N 2^(O(log* N))) operations.

[fdej]

Minor pedantry: the best proved bound for M(N) is O(N log N 8^(log* N)). You can easily pretend that log* N is constant, though...

[Spiritus]

Out of curiosity, what line of work are you in given that you can seemingly pull this stuff from memory? All I read was "lorem ipsum..." until you got to Big O.

[tomp]

Another similar "algorithm": sleep sort.

Basically, take a list of integers, and for each element n create a thread that sleeps n seconds, after which it appends n to the result list.

Assuming the original list is dense (there are no gaps between integers), we can sort it in O(n)!

:)

[robinhouston]

Even if the computation could be done by magic, simply writing the result to memory would take O(n) operations, since the nth Fibonacci number has O(n) bits.

[dastbe]

An issue here is that the the nth Fibonacci number has O(N) bits, so once you move into arbitrary arithmetic calculations your +'s and *'s will become O(N) and O(N logN) respectively.

[isani]

Factoring into powers of 2 seems to me like an unnecessary complication. It's possible to calculate an arbitrary power in O(log N) time without memoization.

  def __get_matrix_power(self, M, p):
    if p == 1:
      return M
    if p % 2 == 1: # odd power
      return self.__multiply_matrices(M, self.__get_matrix_power(M, p - 1))
    else: # even power
      K = self.__get_matrix_power(M, int(p/2))
      return self.__multiply_matrices(K, K)

[nabb]

Also working on pairs of consecutive fibonacci numbers (f_n, f_(n+1)) instead of the matrix [[f_(n+1) f_n] [f_n f_(n-1)]] makes this much simpler.

  def fib(n):
    def fib2(n):
      # returns (f_n, f_(n+1))
      if n == 0:
        return 0, 1
      if n % 2:
        a, b = fib2(n-1)
        return b, a+b
      a, b = fib2(n//2)
      return b*a + a*(b-a), a*a + b*b
    return fib2(n)[0]

[nealwu]

As mentioned, the actual number has O(N) bits and thus can't actually be computed in O(log N) time. So if you just want an approximation, you may as well use Binet's formula. [1]

However, one thing this algorithm can do that Binet's formula can't is compute the N-th Fibonacci number modulo some constant M in O(log N) time. So you could very efficiently compute, say, the last 12 digits of F_123456789. A lot of programming contests (TopCoder, Google Code Jam, etc.) will ask for the answer modulo some constant to take advantage of this trick.

[1] http://en.wikipedia.org/wiki/Fibonacci_number#Closed-form_ex...

[nhaehnle]

A really neat observation here is that Binet's formula and the matrix exponentiation trick are arguably the same thing.

How so? You can think of the operations of Binet's formula as happening in the field Q(sqrt{5}). Since sqrt{5} is a square root, every element of this field has a unique representation of the form a + b sqrt{5} with a and b rational, and so the field can be said to have dimension 2 over the rationals. When you squint at how the field operations work in this representation, you'll find that you end up doing something that looks very similar to the 2x2 matrix multiplications trick; with a bit of rearrangement, it's the same thing - except that the matrix multiplication trick is usually derived via dynamic programming, and Binet's formula is usually derived via power series.

The fact that things often fit together so nicely is perhaps the most beautiful aspect of mathematics.

[wb14123]

I remember there is an equation that could get the Nth fibonacci number directly. (Maybe in the exercises of Concrete Mathematics by Knuth).

[hacknat]

https://en.wikipedia.org/wiki/Fibonacci_number#Closed-form_e...

If ever there was a way to impress an interviewer it is to bust out this baby when asked to make a fibonacci function.

[NhanH]

It's also in SICP, but depend on what you consider the basic step is, it's also O(logN). As someone else mentions, it requires arbitrary precision calculation, and exponential calculation is O(logN).

[msds]

If you diagonalize the matrix in the problem, the formula pops out pretty quickly.

[reikonomusha]

It requires arbitrary precision arithmetic. It is known as Binet's formula.

[richardkiss]

I wrote about this last December:

http://blog.richardkiss.com/?p=398

The matrix math is easier, and the Python code is about a dozen lines.

[lscharen]

This brings back some memories. I remember talking about the different ways of calculating the Fibonacci numbers in an old comp.lang.java usenet thread about recursion, memoization and dynamic programming from 2007,

http://coding.derkeiler.com/Archive/Java/comp.lang.java.prog...

The thread discusses the exponential, linear and logarithmic algorithms.

[sillysaurus3]

In Python:

    def fib(n):
        import math
        r5 = math.sqrt(5)
        x = pow((1 + r5) / 2, n)
        x -= pow((1 - r5) / 2, n)
        return x / r5

  $ python fib.py 4
  3.0

  $ python fib.py 5
  5.0

  $ python fib.py 50
  12586269025.0

  $ python fib.py 5000
  OverflowError: (34, 'Result too large')

So we add arbitrary precision. Uglier, but supports any n:

    def fib(n):
        from math import sqrt
        from decimal import Decimal
        r5 = Decimal(sqrt(5))
        x = pow((1 + r5) / 2, n)
        x -= pow((1 - r5) / 2, n)
        return x / r5

  $ python fib.py 4
  3.000000000000000242931577139

  $ python fib.py 5
  5.000000000000000607328942851

  $ python fib.py 50000
  1.077773489309106593392984005E+10449

  $ python fib.py 50309230
  8.147201098193506535089177522E+10514006

Or does it?

  $ python fib.py 50309230390
  decimal.Overflow: above Emax

This was fun. Computing fib(50309230390) is left as an exercise for the reader.

Wolfram alpha verifies this is correct for Fib 50,000: http://www.wolframalpha.com/input/?i=Fib+50000

[obfuskater]

returning int(round(x)) could act as a bandaid for the slight error in the results.

[webo]

This shows nothing but one's ability to code-up a math formula. The blog post goes into much detail about how to solve the problem using matrices and memoization.

Edit: Looks like sillysaurus3 edited his comment.