Once upon a time, programming interviews would include a fizz-buzz problem to weed out the folks who couldn’t string together a few lines of code. We can debate when and how such things are appropriate interview questions, but one thing that is always appropriate is to use them as inspiration for practising our own skills.¹

There are various common problems offered in such a vein, including fizz-buzz itself, computing certain prime numbers, and computing Fibonacci number. A few years back, I had a go at writing my own Fibonacci function. When I started researching approaches, I discovered an intriguing bit of matrix math, so I learned something while practicing my skills.²

enter the matrix

One problem with calculating a Fibonacci number is that naïve algorithms require n additions. This is obviously expensive for large values of n. But of course, there are some interesting things we can do to improve on this.

In this solution, we observe that we can express the Fibonacci number F(n) using a 2x2 matrix that is raised to the power of n:

[ 1 1 ] n      [ F(n+1) F(n)   ]
[ 1 0 ]    =   [ F(n)   F(n-1) ]

On the face of it, raising someting to the power of n turns n additions into n multiplications. n multiplications sounds worse than n additions, however there is a trick about raising something to a power that we can exploit. Let’s start by writing some code to multiply matrices:

Multiplying two matrices is a little interesting if you have never seen it before:

[ a b ]       [ e f ]   [ ae + bg  af + bh ]
[ c d ] times [ g h ] = [ ce + dg  cf + dh ]

Our matrices always have diagonal symmetry, so we can simplify the calculation because c is always equal to b:

[ a b ]       [ d e ]   [ ad + be  ae + bf ]
[ b c ] times [ e f ] = [ bd + ce  be + cf ]

Now we are given that we are multiplying two matrices with diagonal symmetry. Will the result have diagonal symmetry? In other words, will ae + bf always be equal to bd + ce? Remember that a = b + c at all times and d = e + f provided that each is a power of [[1,1], [1,0]]. Therefore:

ae + bf = (b + c)e + bf
        = be + ce + bf
bd + ce = b(e + f) + ce
        = be + bf + ce

That simplifies things for us, we can say:

[ a b ]       [ d e ]   [ ad + be  ae + bf ]
[ b c ] times [ e f ] = [ ae + bf  be + cf ]

And thus, we can always work with three elements instead of four. Let’s express this as operations on arrays:

[a, b, c] times [d, e, f] = [ad + be, ae + bf, be + cf]

Which we can code in JavaScript:

let times = (...matrices) =>
  matrices.reduce(
    ([a,b,c], [d,e,f]) => [a*d + b*e, a*e + b*f, b*e + c*f]
  );

times([1,1,0]) // => [1, 1, 0]
times([1,1,0], [1,1,0]) // => [2, 1, 1]
times([1,1,0], [1,1,0], [1,1,0]) // => [3, 2, 1]
times([1,1,0], [1,1,0], [1,1,0], [1,1,0]) // => [5, 3, 2]
times([1,1,0], [1,1,0], [1,1,0], [1,1,0], [1,1,0]) // => [8, 5, 3]

To get exponentiation from multiplication, we could write out a naive implementation that constructs a long array of copies of [1,1,0] and then calls times:

let naive_power = (matrix, n) =>
  times(...new Array(n).fill([1,1,0]));

naive_power([1,1,0], 1) // => [1, 1, 0]
naive_power([1,1,0], 2) // => [2, 1, 1]
naive_power([1,1,0], 3) // => [3, 2, 1]
naive_power([1,1,0], 4) // => [5, 3, 2]
naive_power([1,1,0], 5) // => [8, 5, 3]

Very interesting, and less expensive than multiplying any two arbitrary matrices, but we are still performing n multiplications when we raise a matrix to the nth power. What can we do about that?

exponentiation with matrices

Now let’s make an observation: instead of accumulating a product by iterating over the list, let’s Divide and Conquer. Let’s take the easy case: Don’t you agree that times([1,1,0], [1,1,0], [1,1,0], [1,1,0]) is equal to times(times([1,1,0], [1,1,0]), times([1,1,0], [1,1,0]))? And that this saves us an operation, since times([1,1,0], [1,1,0], [1,1,0], [1,1,0]) is implemented as:

times([1,1,0],
  times([1,1,0],
    times([1,1,0],[1,1,0]))

Whereas times(times([1,1,0], [1,1,0]), times([1,1,0], [1,1,0])) can be implemented as:

let double = times([1,1,0], [1,1,0]),
    quadruple = times(double, double);

This only requires two operations rather than three. Furthermore, it is recursive. naive_power([1,1,0], 8) requires seven operations. However, it can be formulated as:

let double = times([1,1,0], [1,1,0]),
    quadruple = times(double, double),
    octuple = times(quadruple, quadruple);

Now we only need three operations compared to seven. Of course, we left out how to deal with odd numbers. Fixing that also fixes how to deal with even numbers that aren’t neat powers of two:

let power = (matrix, n) => {
  if (n === 1) return matrix;

  let halves = power(matrix, Math.floor(n / 2));

  return n % 2 === 0
         ? times(halves, halves)
         : times(halves, halves, matrix);
}

power([1,1,0], 1) // => [1, 1, 0]
power([1,1,0], 2) // => [2, 1, 1]
power([1,1,0], 3) // => [3, 2, 1]
power([1,1,0], 4) // => [5, 3, 2]
power([1,1,0], 5) // => [8, 5, 3]

Now we can perform exponentiation of our matrices, and we take advantage of the symmetry to perform log2n multiplications.

and thus to fibonacci

Thusly, we can write our complete fibonacci function:

let matrixFibonacci = (n) =>
  n < 2
  ? n
  : power([1,1,0], n - 1)[0]

new Array(20).fill(1).map((_, i) => matrixFibonacci(i))
  // => [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181]

We’re done. And this is a win over the typical recursive or even iterative solution for large numbers, because while each operation is more expensive, we only perform log2n operations.³

(This is a translation of a blog post written in 2008)

notes: