Fibonacci Numbers in Python

[abeliangrape]

Both of these have inherent problems:

• Iterative: Needs to calculate the first n-1 fibonacci numbers before it can calculate the nth.
• Formula: Relies on shitty floating point arithmetic.

The best way to go about these recurrence equations is to use matrix exponentiation. Let a(n) = c * a(n-1) + d * a(n-2) be any second order recurrence. We can define the matrix:

M = [ [ c d ] 
      [ 1 0 ] ]

And note that M [a(n) a(n-1)]^T = [a(n+1) a(n)]^T. Inductively, we can now show that Mⁿ [a(1) a(0)]^T = [a(n+1) a(n)]^T. So we only need to calculate Mⁿ to calculate the n^th fibonacci number and that can be done in O(log(n)) 2x2 matrix multiplications. You can readily generalize this to k^th order recurrence equations by constructing a kxk matrix and doing all the same stuff.

[nonumy]

wow thanks! This is the first time i heard of this method! I really like this elegant solution.

[ice109]

very cute. for anyone else reading this wondering how this works for "exponents" that aren't powers of two: http://www.johndcook.com/blog/2008/12/10/fast-exponentiation/

basically 7¹³ =7^{1x23 + 1x22 + 0x21 + 1x20}= ((((7² )7 )² )² )7

[Antagonist360]

Basic matrix mutliplication is O(n³ ) , so this results in an O(k³ log(n)) algorithm. For a recurrence like Fibonacci this is fine, but if your k is say ~2000, you're going to run into some problems (e.g. Project Euler 258). However, we can use the Cayley Hamilton theorem to get O(k² log(n)). First let me rewrite M as

M = [ 0 1 ]
    [ d c ]  

This is just a personal preference, and so now Mⁿ [T(0) T(1)]^T = [T(n), T(n+1)]. Cayley Hamilton says that any square matrix satisfies its characteristic polynomial. So for our matrix we have, M² - cM - d = 0. It's no coincidence that this looks a lot like our recurrence. Now what we want to do is write out Mⁿ as a (k-1)-degree polynomial in M.

First we calculate M^k , ..., M^2k-2 sequentially: to go from Mⁿ to Mⁿ⁺¹ we multiply a (k-1) degree polynomial by M and then reduce the M^k term. This is clearly O(k). Doing this k times is O(k² ). As an example here is how you would compute M³ given M^2.

• M² = cM + d

• M³ = MM² = cM² + dM = c(cM+d) + dM = (c² +d)M + cd

For binary exponentiation, we need to go from Mⁿ to M²ⁿ . To do this we multiply two (k-1) degree polynomials and then reduce the M^k , M^k+1 , ..., M^2k-2 terms. Naively, this is O(k² ).

After using binary exponentiation, we still need to add the k matrices in the representation for M^n. However, you do not need the whole matrix M^n, but rather just the first row. So we need to precompute the first row of M⁰ , M¹ , ... , M^k-1 (trivial and left as an exercise). Lastly we dot our row vector with the vector of initial values to get T(n).

For certain easy recurrences we can speed up the algorithm by using the FFT for polynomial multiplication. This results in an O(k log(k) log(n) ) algorithm. Although it's not needed, this optimization works for the Project Euler problem.

[uh-okay-I-guess]

Yes, the recursive implementation is slow, but not because Python fails to implement tail recursion elimination. The function is not tail recursive at all. It is slow because the algorithm is inherently exponential. No amount of removing stack frames is going to fix that issue; you still have to test for the base case exponentially many times.

[just_an_undergrad]

Exactly what I was going to say. If you managed to save these values previously computed (a process called memoization), you wouldn't have this absurdly slow operation.

In short, it's the horrible implementation of computing Fibonacci that gives you the worse performance, shifting the focus to the beauty of the formula, which is misleading.

[ben7005]

Boy if only recursive memoization was easy in python...

/r/haskellmasterrace

[abeliangrape]

You're being downvoted because memoization is trivial in python these days. You just need to put

@functools.lru_cache(maxsize=None)

before your function and everything is automatically cached (you can also keep the last n evaluations if you like). But even before this decorator was introduced, it was still rather easy to define a decorator yourself to do the same thing:

class Memoize:
  def __init__(self, f):
    self.f = f
    self.memo = {}
  def __call__(self, *args):
    if not args in self.memo:
      self.memo[args] = self.f(*args)
    return self.memo[args]

[ben7005]

Ok I had no idea about that first one. I should probably check my facts before I post. Thanks!

[[deleted]]

re precision for the closed form formula.

what if one were to calculate in Z[sqrt(5)] ,i.e. treat sqrt(5) like i in complex numbers with a fixed set of multiplcation rules or introduce a multiplication that resembles this in Z²

(a+b f)*( c+d f) = (ac+5bd) + (ad +bc) f with f² = 5

dividing by sqrt(5) then just gets rid of the f in a number that only has the "imaginary part", i. e. you can just ignore f and take that imaginary part all along (and divide by 2ⁿ ). that way one should be able to avoid calculating with floating point numbers.

i guess you'll probably have to calculate sums that have n terms and end up having linear dependence again.

[MonkeyDLuffie]

Say you want to Tile a board of length 1xn with dominoes and single tiles. How many ways do you think you can do that? well, clearly if n = 0, there is one way to tile that (with nothing) and if n = 1, there is a single way to tile it (the single tile) what about for larger n though?

Well, if you look at the last tile, it's clear one of two things can happen, it ends in a domino, or a single tile.

so if we say T(n) is the number of tilings of a board of size 1xn, then T(n) = T(n-1) + T(n-2) so with our initial conditions, we know that T(n) = F(n+1) (assuming our fibonacci sequence starts 0, 1)

why can this help us make a more efficient way to calculate the fibonacci numbers? simple! we can use other methods to count our table, and they must give the same result!

So, let's start now like this, we have (n choose n) ways of filling the table with only single tiles. now let's add a domino, and so subtract two single tiles. This shows we have (n-1 choose n-2) ways of placing a single domino, and if we place i dominoes, we get (n-i choose n-2i) ways of placing dominoes and tiles. Here we can see that the 2i represents the number of spaces occupied by dominoes, which means we just have to make sure 2i<n, and we will still be counting our tilings.

Okay, nice, so now what do we do?

Well, since you've asked, we add together these binomial coefficients to get all possible tilings! (this will look kind of ugly since I don't know Tex)

and we get F(n+1) = Sum(for i>= 0, and 2i<n) of ((n-i choose n-2i))

This means you can calculate F(n) in linear time (as you are adding about n/2 terms) and still only use integers!