Calculating Fibonacci Numbers, Quickly and Exactly

Richard

December 28, 2013

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 \\ F_{2n+1} & = F_{n}^2 + F_{n+1}^2 \end{align} $$

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.