# 0 1 1 2 3 5 8 13…

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 Fn = (φn-ψn)/√5 Ln = φn+ψn```

We can now derive:

```5Fn2 = φ2n+ψ2n-2(-1)n = L2n-2(-1)n Ln2 = φ2n+ψ2n+2(-1)n = L2n+2(-1)n LnFn = (φ2n-ψ2n)/5 = F2n```

We can also use the formulas to prove:

`2φn = Ln + Fn√5`

which gives us:

`Ln+m + Fn+m√5 = 2φn+m = 2φnφm = (LnLm+5FnFm)/2 + (LnFm+FnLm)√/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`):

```Ln+m = (LnLm+5FnFm)/2 Fn+m = (LnFm+FnLm)/2```

Two special cases for this:

```L2n = (Ln2+5Fn2)/2 F2n = LnFn```

```Ln+1 = (Ln+5Fn)/2 Fn+1 = (Ln+Fn)/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 Ln 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:

`Fn+m = Fn+1Fm+FnFm-1`

which gives us:

```F2n = FnFn+2FnFn-1 F2n-1 = FnFn+Fn-1Fn-1```

and

```F2n+1 = Fn+1Fn+1+FnFn F2n = 2FnFn+1-FnFn```

```
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 `(an+2,an+1)` is `M(an+1,an)` for 2×2 matrix `M = (1 1 1 0)`. We then have `(Fn+1,Fn) = Mn(1 0)` and we can use conventional matrix operations, including inversion, to compute `Fn`. We can apply this to any linear recurrence, eg. `an+2 = Aan+1+Ban` corresponds to the matrix `(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 `(Ln+1,Ln) = Mn(1 2)`, 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 `(a+b a a b)` (cf. HAKMEM: http://www.inwap.com/pdp10/hbaker/hakmem/recurrence.html). But that probably needs another post.