the best fibonacci implementation?

I’ve been meaning to write about this ever since I came across it in Abelson & Sussman’s wizard book.

I mean, Fibonacci is the classic example of a recursive function. And it’s often just implemented as:

def fib(n):
    if n == 0 or n == 1:
        return 1
    return fib(n-1) + fib(n-2)
print(fib(5))

8

And while this seems all well and good on the surface, it turns out that this is far from an efficient implementation of the sequence. In particular, we end up computing calls to fib that we’ve already computed. In other words, we perform a lot of redundant computation.

For example, in evaluating fib(5) we end up calling fib(4) and fib(3). But notice how in evaluating fib(4) we call fib(3) again, alongside fib(2). This problem gets worse for larger values of n.

I would show you in Python, but it’s not really built for recursion and you’ll quickly exceed the maximum recursion depth.

You could technically circumvent this, but why bother.

Here’s the same implementation of our original Python fib function in Scheme:

(define (fib n)
  (if (or (= n 0) (= n 1))
      1
      (+ (fib (- n 1)) (fib (- n 2)))))
(fib 5)

8

Prefix notation aside, to get a sense of how bad this implementation of Fibonacci is in practice, even in a language that takes recursion seriously, it takes just under 7s to compute (fib 35) on my machine ¹.

That’s really not great.

Fortunately, we can do much better. And we don’t even need to redefine the way we compute the \(n^{\text{th}}\) Fibonacci to do it!

Thankfully, it turns out that there is a clever way of mitigating redundant computations in recursive functions: memoisation.

In essence, every time we make a call to a recursive function we store its result in a table indexed by the arguments of that function call. Put simply, we keep a record of which arguments give us which results when we call our recursive function.

We do this so that each time we attempt a call to our recursive function we check this record to see if we’ve already previously called our function with this set of arguments. If we have, we just return the result already stored in the record, the one we put there the first time we evaluated that function with those arguments. If not, we evaluate the function call and record the result.

It sounds abstract but in our case it’s really just about preventing us from shooting ourselves in the foot by having to evaluate (fib 3) from scratch, when we end up calling it in the evaluation of (fib 5) and (fib 4), for example.

Anyways…

Once we memoise our procedure² we can compute (fib 35) in so little time it might as well be 0s.

It takes .01s to evaluate (fib 1000) – which is 209 digits long by the way.³ (We can even compute (fib 10000) and its 2090 digits in ~1.33s)

But that’s not cause for celebration – we can do better!

How much better you may ask?

Well, if you care at all about time complexity, know that our first implementation was in exponential time, this memoised version is in linear time, and we can get it all the way down to logarithmic time!

Better yet, we can get the space complexity down to constant space!

(This last part is important because if you wanted to compute (fib 100000), for example. Even with our memoised version you would probably run out of memory, at least I did.)

So here it is – the best Fibonacci implementation! (maybe?)

cleverness + tail recursion

Here’s the thing, computing Fibonacci more efficiently requires us to change the way we think about the sequence in the first place.

More specifically, we need to abandon our traditional understanding of “the \(n^{\text{th}}\) Fibonacci number is the sum of the \((n-1)^{\text{th}}\) and \((n-2)^{\text{th})}\) (nd?)”, in favour of something a little cleverer.

Look at the Fibonacci sequence:

0 1 1 2 3 5 8 13…

and consider it as successive, overlapping pairs:

(0,1) (1,1) (1,2) (2,3) (3,5) (5,8) (8,13) …

Doing so may seem odd at first, and in some ways it is. But the point is that an interesting pattern emerges when you consider the relationship between these specific pairs (i.e. how the elements in the pairs change from pair to pair).

The easiest thing to notice is that the second element of a pair becomes the first element of the pair that follows it. This is a direct consequence of our pairs being overlapping. In other words, this particular pattern is given by design.

We can formalise this observation by saying that given a pair \((a,b)\) in our sequence, we know the pair that follows it will look like \((b,?)\).

The question is then whether there’s some kind of relationship between a pair and the second element of the pair that follows it?

Well, take a look at just the second elements of each pair in our sequence:

1 1 2 3 5 8 13…

…the Fibonacci sequence!

Again, this is by design, by definition of our sequence of pairs.

Look:

• (0,1) is the first and second element of the Fibonacci sequence
• (1,1) is the second and third element of the Fibonacci sequence
• (1,2) is the third and fourth element of the Fibonacci sequence
• (2,3) is the fourth and fifth element of the Fibonacci sequence
• …

What does this mean though?

This means that the second element of a pair is the sum of the two previous Fibonacci numbers, but this is just the previous pair in the sequence of pairs!

Using our previous formalisation this means that a given pair \((a,b)\) will be followed by the pair \((b,a+b)\).

Another way to express this is that the elements \(a\) and \(b\) in a pair \((a,b)\) get transformed into the elements \(b\) and \(a+b\) respectively. We can use the notation \(a \leftarrow b\) and \(b \leftarrow a + b\) to express this.

We can even just implement Fibonacci with this:

def fib(n):
    a = 0
    b = 1
    while n > 0:
        temp_a = a
        a += b
        b = temp_a
        n -= 1
    return a
for i in range(8):
    print(fib(i))

0
1
1
2
3
5
8
13

But that’s not it. This is not “the best Fibonacci” implementation, not yet.

It turns out that the transformation that describes \(a \leftarrow b\) and \(b \leftarrow a + b\) can be understood as a specific instance of a more general transformation of the form \(a \leftarrow ap + bq\) and \(b \leftarrow aq + bq + bp\) where \(q=1\) and \(p=0\) (just substitute those values to check we get back our original transformation).

For the sake of convenience, we can call this more general transformation \(T_{pq}\).

And the reason we do all this is because it allows us to ask whether applying the \(T_{pq}\) transformation twice is itself a kind of transformation \(T_{p’q'}\), where \(p'\) and \(q'\) are perhaps different to \(p\) and \(q\)?

While you can probably answer this question through repeated substitution between \(a \leftarrow ap + bq\) and \(b \leftarrow aq + bq + bp\), I personally find this approach easily confusing.

Instead, I prefer representing \(T_{pq}\) as a matrix such that

\begin{equation*} T_{pq}\begin{pmatrix}a \\ b \end{pmatrix} = \begin{pmatrix} ap + bq \\ aq + bq + bp \end{pmatrix}. \end{equation*}

And if you know anything about matrix operations, you can figure out \(T_{pq}\)’s matrix by inspection. If not just know that it happens to be \(\begin{pmatrix} p & q \\ q & p+q \end{pmatrix}\).

Given this, linear algebra tells us that applying \(T_{pq}\) twice is just squaring it:

\begin{align*} T_{pq}^{2} &= \begin{pmatrix} p & q \\ q & p+q \end{pmatrix} \begin{pmatrix} p & q \\ q & p+q \end{pmatrix}\\ &= \begin{pmatrix} p^{2}+q^{2} & pq + q(p+q) \\ pq + q(p+q) & q^{2} + (p+q)^{2} \end{pmatrix} \\ &= \begin{pmatrix} p^{2}+q^{2} & pq + pq+q^{2} \\ pq + pq+q^{2} & q^{2} + p^{2}+2pq+q^{2} \end{pmatrix} \\ &= \begin{pmatrix} p^{2}+q^{2} & 2pq+q^{2} \\ 2pq +q^{2} & 2q^{2} + 2pq + p^{2}\end{pmatrix} \end{align*}

Ok. But does this look like another kind of \(T_{pq}\) with different/new values for \(p\) and \(q\), i.e. \(p'\) and \(q'\).

Again, some careful inspection tells us that \(p'=p^{2}+q^{2}\) and \(q'=2pq+q^{2}\).

We’ve done a fair at this point. So let’s recap:

• We can represent the Fibonacci sequence as a somewhat odd sequence of overlapping pairs.
• Turns out that there is a transformation that describes how to turn any one of these pairs in into the pair that follows it.
• That transformation can be understood as an instance of a more general transformation \(T_{pq}\).
• Applying that transformation twice is equivalent to another instance of the general transformation \(T_{pq}\), i.e. \(T_{p’q'}\).
• Hence, applying \(T_{p’q'}\) to a pair in our sequence won’t give us back the pair that follows it, that’s \(T_{pq}\), instead, it will give us the pair that comes after that. In other words, \(T_{p’q'}\) “skips” pairs in our sequence, generating every other pair.

Ok. But why all this work?

Because, we can now reframe our problem of computing the \(n^{\text{th}}\) Fibonacci as applying \(n\) transformations \(T_{pq}\) (i.e. \(T_{pq}^n\)) to the original pair (0,1) of our sequence. And knowing \(T_{pq}^{2}\) allows us to leverage, the finale piece of this puzzle:

successive squaring

The best way to understand successive squaring is as a clever way of quickly calculating some base \(b\) to the power \(n\), i.e. \(b^{n}\). Instead of just multiplying \(b\) by itself \(n\) times, performing \(n\) multiplications, you can successively square \(b\): \(b,b^{2},b^{4}\ldots\), getting you to \(b^{n}\) in much fewer steps (\(\log_2 n\) steps to be more precise).

It’s all thanks to the observation that \(b^n\) can be recursively defined as:

• \(b^{n} = (b^{n/2})^{2}\) if \(n\) is even
• \(b^{n} = b \cdot b^{n-1}\) if \(n\) is odd

(define (exp b n)
  (cond ((= n 0) 1)
        ((= n 1) b)
        ((= (remainder n 2) 0) (square (exp b (/ n 2))))
        (else (* b (exp b (- n 1))))))
(exp 3 123)

;Value: 48519278097689642681155855396759336072749841943521979872827

That’s it!

The only difference in our case is that we are looking to raise a transformation, and not a number, to the power of \(n\), i.e \(b= T_{pq}\). But this makes no difference since we’ve already figured out how to square our \(T_{pq}\) transformation – you transform the transformation by \(T_{p’q'}\) (i.e. you turn \(p\) into \(p'\) and \(q\) into \(q'\))!

All that’s left is to implement the final solution and bask in its performance:

(define (even? n) (= (remainder n 2) 0))
; calculates the transformation and its application to a starting pair (a,b)
(define (exp-trans a b p q n)
  (cond ((= n 0) a)
        ((even? n)
         (exp-trans a
                    b
                    (+ (square p) (square q)) ;; T_p'q'
                    (+ (* 2 p q) (square q)) ;; T_p'q'
                    (/ n 2)))
        (else (exp-trans (+ (* a p) (* b q)) ;; T_pq
                         (+ (* a q) (* b q) (* b p)) ;; T_pq
                         p
                         q
                         (- n 1)))))
(define (fib n) (exp-trans 0 1 0 1 n))
(define (fib-digits n) (string-length (number->string (fib n))))
(time fib-digits 1000000)

Elapsed run time: .72s
Garbage collector run time: 0.s
Elapsed real time: .72s
;Value: 208988

PS: what about tail recursion?

It’s true, we did name drop tail recursion in one of the headings without really explaining it, nor how it contributes to the “best” Fibonacci implementation we arrived at.

This is in large part because it only influences the space complexity of our implementation, not its time complexity (our main focus). As in, it doesn’t really have anything to do with the speed of the algorithm we developed, it’s more of a nice language feature that allows our algorithm to consume very little space when running.

So while there wasn’t really a good place to go on a tangent like tail recursion, you can read more about it in Section 1.2.1 of the wizard book.