That formula fails for fib(71) with double precision floating point arithmetic:

    In [1]: def fib(n):
       ...:     a,b = 0,1
       ...:     for i in range(0,n): a,b = b,a+b
       ...:     return a
       ...: 

    In [2]: def fib2(n):
       ...:     phi = (1+sqrt(5))/2
       ...:     return round(phi**n / sqrt(5))
       ...: 

    In [3]: fib2(71)
    Out[3]: 308061521170130.0

    In [4]: fib(71)
    Out[4]: 308061521170129L

And it's not that the latter is not exactly representable as a floating point number:

    In [5]: float(_)
    Out[5]: 308061521170129.0

So you'd have to use bigger precision numbers. How many digits are you going to use? You have to use enough not to produce a wrong result and not too many to slow down the algorithm. I bet that the resulting algorithm is slower than the standard [1 1; 1 0]^n matrix exponentiation. Note that your formula comes from the eigenvalue expansion of this matrix exponentiation, which -- despite what you might hear from undergraduate linear algebra courses -- is not a good method for matrix exponentiation. On the contrary, finding eigenvalues is done by a process analogous to matrix exponentiation: http://en.wikipedia.org/wiki/Power_iteration