The well-known Fibonacci series \(F_n\) can be defined as follows:

\(F_n =

\begin{cases}

0 & n = 0 \\

1 & n = 1 \\

F_{n-2} + F_{n-1} & n \ge 2\\

\end{cases}\)

Let’s use a few facts about matrices to find a quick way to calculate terms in this famous series.

**Lemma**

Let \(A = \begin{bmatrix}

0 & 1 \\

1 & 1

\end{bmatrix}\). Then \(A^n = \begin{bmatrix} F_{n-1} & F_n \\ F_n & F_{n+1} \end{bmatrix}\).

**Proof**

The \(n = 1\) case follows immediately from the definitions of \(A \text{ and } F_n\).

Suppose the statement is true for n. Then

$$\begin{align}

A^{n+1} & = A^n A \\

& = \begin{bmatrix} F_{n-1} & F_n \\ F_n & F_{n+1} \end{bmatrix} \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix} \\

& = \begin{bmatrix} F_n & F_{n-1} + F_n \\ F_{n+1} & F_n + F_{n+1} \end{bmatrix} \\

& = \begin{bmatrix} F_n & F_{n+1} \\ F_{n+1} & F_{n+2} \end{bmatrix}

\end{align}$$

And our result follows by induction. QED

Now, notice that

\(\begin{align}

\begin{bmatrix} F_{2n-1} & F_{2n} \\ F_{2n} & F_{2n+1} \end{bmatrix} & = A^{2n} \\

& = A^n A^n \\

& = \begin{bmatrix} F_{n-1} & F_n \\ F_n & F_{n+1} \end{bmatrix} \begin{bmatrix} F_{n-1} & F_n \\ F_n & F_{n+1} \end{bmatrix} \\

& = \begin{bmatrix} F_{n-1}^2 + F_n ^ 2 & F_{n-1} F_n + F_n F_{n+1} \\ F_{n-1} F_n + F_n F_{n+1} & F_n^2 + F_{n+1} ^ 2 \end{bmatrix}

\end{align}\)

Using this identity, we can write \(F_{2n}\) and \(F_{2n+1}\) in terms of \(F_n\) and \(F_{n+1}\).

\(\begin{align}

F_{2n} & = F_{n-1} F_n + F_n F_{n+1} \\

& = (F_{n+1} – F_n) F_n + F_n F_{n+1} \\

& = F_n F_{n+1} – F_n^2 + F_n F_{n+1} \\

& = 2 F_{n+1} F_n – F_n ^ 2

\end{align} \)

\(F_{2n+1} = F_{n}^2 + F_{n+1}^2\)

We now have a way of calculating \(F_{2n}\) and \(F_{2n+1}\) by calculating only a few of the smaller terms in the sequence. This yields some wildly efficient Python code:

```
```def fib2(N):
if N == 0: return (0, 1)
half_N, is_N_odd = divmod(N, 2)
f_n, f_n_plus_1 = fib2(half_N)
f_n_squared = f_n * f_n
f_n_plus_1_squared = f_n_plus_1 * f_n_plus_1
f_2n = 2 * f_n * f_n_plus_1 - f_n_squared
f_2n_plus_1 = f_n_squared + f_n_plus_1_squared
if is_N_odd:
return (f_2n_plus_1, f_2n + f_2n_plus_1)
return (f_2n, f_2n_plus_1)
def fib(N):
return fib2(N)[0]

And it’s fast! It’s \(O(\log N)\) in fact:

$ pypy -m timeit -s 'import expfib' 'expfib.fib(500000)' 10 loops, best of 3: 22.4 msec per loop

Compare to the naïve, iterative version, which is \(O(N)\):

```
```def fib2_linear(N):
f_n, f_n_plus_1 = (0, 1)
while N > 0:
N -= 1
f_n, f_n_plus_1 = f_n_plus_1, f_n + f_n_plus_1
return (f_n, f_n_plus_1)
def fib_linear(N):
return fib2_linear(N)[0]

$ pypy -m timeit -s 'import myfib' 'myfib.fib_linear(500000)' 10 loops, best of 3: 2.76 sec per loop

This post was inspired by a post by Lee Phillips and as an excuse to play around with MathJax. The best guide I found for getting started is on Stack Exchange.

Great code. I love it!

Neat.

It took me a little effort to convince myself that your final algorithm is logarithmic and not linear in time. A naïve implementation would be linear, but in a more interesting way than usual: the call tree would have a branch factor of two but would now only be log2(n) deep.

But I see that you take advantage of the two right-hand-side terms always being adjacent to reduce the branch factor to one. Very nice.