For some reason I spent some time over Christmas looking at algorithms for computing Fibonacci numbers, so here’s the first post for 2013.

The Binet formulas for the Fibonacci and Lucas numbers are:

`φ = (1+√5)/2`

ψ = (1-√5)/2

φψ = -1

F_{n} = (φ^{n}-ψ^{n})/√5

L_{n} = φ^{n}+ψ^{n}

We can now derive:

`5F`

_{n}^{2} = φ^{2n}+ψ^{2n}-2(-1)^{n} = L_{2n}-2(-1)^{n}

L_{n}^{2} = φ^{2n}+ψ^{2n}+2(-1)^{n} = L_{2n}+2(-1)^{n}

L_{n}F_{n} = (φ^{2n}-ψ^{2n})/5 = F_{2n}

We can also use the formulas to prove:

`2φ`

^{n} = L_{n} + F_{n}√5

which gives us:

`L`

_{n+m} + F_{n+m}√5 = 2φ^{n+m} = 2φ^{n}φ^{m} = (L_{n}L_{m}+5F_{n}F_{m})/2 + (L_{n}F_{m}+F_{n}L_{m})√/2

So, equating the rational and irrational parts (as we may – if `a+bk = c+dk`

with `a`

,`b`

,`c`

,`d`

rational, then `k`

is rational unless `a=c`

and `b=d`

):

`L`

_{n+m} = (L_{n}L_{m}+5F_{n}F_{m})/2

F_{n+m} = (L_{n}F_{m}+F_{n}L_{m})/2

Two special cases for this:

`L`

_{2n} = (L_{n}^{2}+5F_{n}^{2})/2

F_{2n} = L_{n}F_{n}

`L`

_{n+1} = (L_{n}+5F_{n})/2

F_{n+1} = (L_{n}+F_{n})/2

With these recurrences, we can write some code:

def fib1aux(n): """Return F(n),L(n)""" if n == 0: return 0,2 else: a,b = fib1aux(n//2) a,b = a*b,(5*a*a+b*b)//2 if n%2 == 0: return a,b else: return (a+b)//2, (5*a+b)//2 return a,b def fib1(n): return fib1aux(n)[0]

We can save a bignum multiplication in each recursive step by replacing `a,b = a*b,(5*a*a+b*b)//2`

with something like `ab = a*b; t = (a+b)*(5*a+b); a,b = ab,t//2-3*ab`

(suggestion due to Robin Houston at http://bosker.wordpress.com/2011/07/27/computing-fibonacci-numbers-using-binet%E2%80%99s-formula/).

Alternatively, we can follow the suggestion of Knuth in TAOCP (solution to exercise 4.6.3.26) and define:

def fib2aux(n): """Return F(n),L(n)""" if n == 0: return 0,2 else: a,b = fib2aux(n//2) a,b = a*b,n//2%2*4-2+b*b if n%2 == 0: return a,b else: return (a+b)//2, (5*a+b)//2 def fib2(n): return fib2aux(n)[0]

In both cases, we benefit from avoiding the final calculation of L_{n} and define something like:

def fib2(n): if (n%2 == 0): a,b = fib2aux(n//2) return a*b else: a,b = fib2aux(n-1) return (a+b)//2

or, for a non-recursive solution (stealing some ideas from Robin Houston again):

def bits(n): r = [] while n!=0: r.append(n%2) n >>=1 r.reverse() return r def fib3(n): a,b = 0,2 last = 0 for i in bits(n//2): a,b = a*b,last*4-2+b*b last = i if i: a,b = (a+b)//2, (5*a+b)//2 if n%2 == 0: return a*b else: a,b = a*b,last*4-2+b*b return (a+b)//2

Of course, we don’t have to use Lucas numbers, a pleasant recurrence for Fibonacci numbers alone is:

`F`

_{n+m} = F_{n+1}F_{m}+F_{n}F_{m-1}

which gives us:

`F`

_{2n} = F_{n}F_{n}+2F_{n}F_{n-1}

F_{2n-1} = F_{n}F_{n}+F_{n-1}F_{n-1}

and

`F`

_{2n+1} = F_{n+1}F_{n+1}+F_{n}F_{n}

F_{2n} = 2F_{n}F_{n+1}-F_{n}F_{n}

leading to:

def fibaux4(n): if n == 0: return 0,1 elif n == 1: return 1,0 else: a,b = fibaux4((n+1)//2) if (n%2 == 0): return a*(a+2*b),a*a+b*b else: return a*a+b*b,b*(2*a-b) def fib4(n): return fibaux4(n)[0]

with an elegant symmetry between the odd and even cases. This seems to be solution suggested by Dijkstra in http://www.cs.utexas.edu/~EWD/ewd06xx/EWD654.PDF.

There is, of course more: notably the matrix based solutions – we observe that `(a`

is _{n+2},a_{n+1)}`M(a`

for 2×2 matrix _{n+1},a_{n})`M = (1 1 1 0)`

. We then have `(F`

and we can use conventional matrix operations, including inversion, to compute _{n+1},F_{n}) = M^{n}(1 0)`F`

. We can apply this to any linear recurrence, eg. _{n}`a`

corresponds to the matrix _{n+2} = Aa_{n+1}+Ba_{n}`(A B 1 0)`

.

It’s notable that `M`

only depends on the recurrence and not on the initial values, so we can compute, for example, the Lucas numbers with `(L`

, and there is a regularity to the structure of M that can be used to optimise the computation, eg. for Fibonacci, M is of the form _{n+1},L_{n}) = M^{n}(1 2)`(a+b a a b)`

(cf. HAKMEM: http://www.inwap.com/pdp10/hbaker/hakmem/recurrence.html). But that probably needs another post.