Benchmark Fibonacci

fib.py

from timeit import timeit

class MatrixFibonacci:
    Q = [[1, 1],
         [1, 0]]

    def __init__(self):
        self.__memo = {}

    def __multiply_matrices(self, M1, M2):
        """Matrices miltiplication
        (the matrices are expected in the form of a list of 2x2 size)."""

        a11 = M1[0][0]*M2[0][0] + M1[0][1]*M2[1][0]
        a12 = M1[0][0]*M2[0][1] + M1[0][1]*M2[1][1]
        a21 = M1[1][0]*M2[0][0] + M1[1][1]*M2[1][0]
        a22 = M1[1][0]*M2[0][1] + M1[1][1]*M2[1][1]
        r = [[a11, a12], [a21, a22]]
        return r

    def __get_matrix_power(self, M, p):
        """Matrix exponentiation (it is expected that p that is equal to the power of 2)."""

        if p == 1:
            return M
        if p in self.__memo:
            return self.__memo[p]
        K = self.__get_matrix_power(M, int(p/2))
        R = self.__multiply_matrices(K, K)
        self.__memo[p] = R
        return R

    def get_number(self, n):
        """Getting the nth Fibonacci number
        (a non-negative integer number is expected as n)."""
        if n == 0:
            return 0
        if n == 1:
            return 1
        # Factoring down the passed power into the powers that are equal to the power of 2), 
        # i.e. 62 = 2^5 + 2^4 + 2^3 + 2^2 + 2^0 = 32 + 16 + 8 + 4 + 1.
        powers = [int(pow(2, b))
                  for (b, d) in enumerate(reversed(bin(n-1)[2:])) if d == '1']
        # The same, but less pythonic: http://pastebin.com/h8cKDkHX
        
        matrices = [self.__get_matrix_power(MatrixFibonacci.Q, p)
                    for p in powers]
        while len(matrices) > 1:
            M1 = matrices.pop()
            M2 = matrices.pop()
            R = self.__multiply_matrices(M1, M2)
            matrices.append(R)
        return matrices[0][0][0]

def fib1(n):
    a, b = 0, 1
    for _ in xrange(n):
        a, b = b, a+b
    return a

def fib2(n):
    matrix = MatrixFibonacci()
    return matrix.get_number(n)

def fib3(n):
    def mul(A, B):
        return (A[0]*B[0] + A[1]*B[2], A[0]*B[1] + A[1]*B[3], 
                A[2]*B[0] + A[3]*B[2], A[2]*B[1] + A[3]*B[3])

    def powm(X, n):
        if n==0: return (1, 0, 0, 1)
        if n==1: return X
        A = powm(X, n>>1)
        A = mul(A, A)
        if n&1: A = mul(A, X)
        return A  

    return powm((1, 1, 1, 0), n)[1]

for i in xrange(0, 500000, 10000):
    print i, '\t', timeit(lambda: fib1(i), number=1), '\t', timeit(lambda: fib2(i), number=1), '\t', timeit(lambda: fib3(i), number=1)