Calculating nth fibonacci number using the formulae in python

Question

I am calculating the n-th fibonacci number using (a) a linear approach, and (b) this expression

Python code:

'Different implementations for computing the n-th fibonacci number'

def lfib(n):
    'Find the n-th fibonacci number iteratively'
    a, b = 0, 1
    for i in range(n):
        a, b = b, a + b
    return a

def efib(n):
    'Compute the n-th fibonacci number using the formulae'
    from math import sqrt, floor
    x = (1 + sqrt(5))/2
    return long(floor((x**n)/sqrt(5) + 0.5))


if __name__ == '__main__':
  for i in range(60,80):
    if lfib(i) != efib(i):
      print i, "lfib:", lfib(i)
      print "   efib:", efib(i)

For n > 71 I see that the two functions return different values.

Is this due to floating point arithmetic involved in efib()? If so, is it then advisable to calculate the number using the matrix form?

Answer 1

You are indeed seeing rounding errors.

The matrix form is the more accurate and much faster algorithm. Literateprograms.org lists a good implementation, but it also lists the following algorithm based on Lucas numbers:

def powLF(n):
    if n == 1:     return (1, 1)
    L, F = powLF(n//2)
    L, F = (L**2 + 5*F**2) >> 1, L*F
    if n & 1:
        return ((L + 5*F)>>1, (L + F) >>1)
    else:
        return (L, F)

def fib(n):
    if n & 1:
        return powLF(n)[1]
    else:
        L, F = powLF(n // 2)
        return L * F

Take a look at Lecture 3 of the MIT Open Courseware course on algorithms for a good analysis of the matrix approach.

Both the above algorithm and the matrix approach has Θ(lg n) complexity, just like the naive recursive squaring method you used, yet without the rounding problems. The Lucas numbers approach has the lowest constant cost, making it the faster algorithm (about twice as fast as the matrix approach):

>>> timeit.timeit('fib(1000)', 'from __main__ import fibM as fib', number=10000)
0.40711593627929688
>>> timeit.timeit('fib(1000)', 'from __main__ import fibL as fib', number=10000)
0.20211100578308105

Answer 2

Is this due to floating point arithmetic involved in efib()?

Yes, it is. Within efib you have

>>> log(x**72)/log(2)
49.98541778140445

and Python floats have about 53 bits of precision on x86-64 hardware, so you're running close to the edge.