lin-alg-fib.py

import sys
import numpy as np
import operator

# def discriminant(a,b,c):
    # return (b**2) - 4*a*c

# def quadratic(a,b,c):
    # """Compute roots from a polynomial of the form ax^2 + bx + c."""
    # ops = [operator.add, operator.sub]
    # return [op(-b, discriminant(a,b,c)**(.5))/(2*a) for op in ops]

def unwrap_R(R):
    """Going to look like:

    [[ fib(n-1), fib(n) ]
     [ fib(n),   fib(n) ]]

    so we can take any of the values except for (0,0).
    """
    return int(R[1,1])

def fib(n):
    """Compute fibonacci number using linear algebra."""
    
    # Magic matrix that represents the fibonacci function
    R = np.array([[0,1],[1,1]])

    # Constants computed from quadratic formula x^2 - x - 1 = 0
    phi = .5 - 5**(.5)/2
    phi_prime = .5 + 5**(.5)/2

    # Or:
    # phi, phi_prime = quadratic(1., -1., -1.)

    # Formula derived at end of Eugene's talk
    prefix = 1 / (phi_prime - phi)
    expr = (phi_prime**n - phi**n)*R + (phi_prime**(n-1) - phi**(n-1))
    return unwrap_R(prefix * expr)

if __name__ == '__main__':
    print(fib(int(sys.argv[1])))