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.