matthew arcus

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 = (φnn)/√5
Ln = φnn

We can now derive:

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

We can also use the formulas to prove:

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
        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

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

def fib2aux(n):
    """Return F(n),L(n)"""
    if n == 0: return 0,2
        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
        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:
        n >>=1
    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
        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


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

leading to:

def fibaux4(n):
    if n == 0: return 0,1
    elif n == 1: return 1,0
        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‎.

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: But that probably needs another post.