In [1]:

using LinearAlgebra

Fibonacci recurrence¶

The Fibonacci numbers are:

$$ 1,1,2,3,5,8,13,21,34,\ldots $$

Each number $f_n$ in the sequence is the sum of the previous two, defining the recurrence relation:

$$ f_n = f_{n-1} + f_{n-2} $$

Perhaps the most obvious way to implement this in a programming language is via recursion:

In [2]:

function slowfib(n)
    if n < 2
        return BigInt(1) # use bigint type to support huge integers
    else
        return slowfib(n-1) + slowfib(n-2)
    end
end

Out[2]:

slowfib (generic function with 1 method)

Note that there is a slight catch: we have to make sure to do our computations with the BigInt integer type, which implements arbitrary precision arithmetic. The Fibonacci numbers quickly get so big that they overflow the maximum representable integer using the default (fast, fixed numbrer of binary digits) hardware integer type.

In [3]:

[slowfib(n) for n = 1:10]

Out[3]:

10-element Vector{BigInt}:
  1
  2
  3
  5
  8
 13
 21
 34
 55
 89

Not that it matters for toy calculations like this, but there are much faster ways to compute Fibonacci numbers than the recursive function defined above. The GMP library used internally by Julia for BigInt arithmetic actually provides an optimized Fibonacci-calculating function mpz_fib_ui that we can call if we want to using the low-level ccall technique:

In [4]:

function fastfib(n)
    z = BigInt()
    ccall((:__gmpz_fib_ui, :libgmp), Cvoid, (Ref{BigInt}, Culong), z, n)
    return z
end

Out[4]:

fastfib (generic function with 1 method)

In [5]:

[fastfib(i) for i = 1:100]

Out[5]:

100-element Vector{BigInt}:
                     1
                     1
                     2
                     3
                     5
                     8
                    13
                    21
                    34
                    55
                    89
                   144
                   233
                     ⋮
   1779979416004714189
   2880067194370816120
   4660046610375530309
   7540113804746346429
  12200160415121876738
  19740274219868223167
  31940434634990099905
  51680708854858323072
  83621143489848422977
 135301852344706746049
 218922995834555169026
 354224848179261915075

It's about 1000x faster even for the 20th Fibonacci number. It turns out that the recursive algorithm is pretty terrible — the time increases exponentially with n.

In [6]:

@time fastfib(20)
@time fastfib(20)
@time fastfib(20)
@time slowfib(20)
@time slowfib(20)
@time slowfib(20)

  0.000007 seconds (2 allocations: 40 bytes)
  0.000002 seconds (2 allocations: 40 bytes)
  0.000002 seconds (2 allocations: 40 bytes)
  0.002284 seconds (43.78 k allocations: 940.625 KiB)
  0.002528 seconds (43.78 k allocations: 940.625 KiB)
  0.002812 seconds (43.78 k allocations: 940.625 KiB)

Out[6]:

Fibonacci as matrices¶

We can represent the Fibonacci recurrence as repeated multiplication by a $2 \times 2$ matrix, since:

$$ \begin{pmatrix} f_{n+1} \\ f_n \end{pmatrix} = \underbrace{\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}}_F \begin{pmatrix} f_{n} \\ f_{n-1} \end{pmatrix} $$

In [7]:

F = [1 1
     1 0]

Out[7]:

2×2 Matrix{Int64}:
 1  1
 1  0

In [8]:

F^7 * [1,1]

Out[8]:

2-element Vector{Int64}:
 34
 21

So, plugging in $f_1 = 1, f_2 = 1$, then

$$ \begin{pmatrix} f_{n+2} \\ f_{n+1} \end{pmatrix} = F^n \begin{pmatrix} 1 \\ 1 \end{pmatrix} $$

and the key to understanding $F^n$ is the eigenvalues of $F$:

In [9]:

eigvals(F)

Out[9]:

2-element Vector{Float64}:
 -0.6180339887498948
  1.618033988749895

Analytically, we can easily solve this $2 \times 2$ eigenproblem to show that the eigenvalues are $(1 \pm \sqrt{5})/2$ (just the roots of the quadratic characteristic polynomial $\det (F-\lambda I) = \lambda^2 - \lambda - 1$):

In [10]:

(1 + √5)/2

Out[10]:

1.618033988749895

In [11]:

(1 - √5)/2

Out[11]:

-0.6180339887498949

For example, to compute $f_{100}$, we can multiply $F^{98}$ by $(1,1)$ (again converting to BigInt using big first to avoid overflow):

In [12]:

big.(F)^98 * [1, 1]

Out[12]:

2-element Vector{BigInt}:
 354224848179261915075
 218922995834555169026

This matches our fastfib function from above:

In [13]:

fastfib(100), fastfib(99)

Out[13]:

(354224848179261915075, 218922995834555169026)

An important thing about $F^n$ is that, for large $n$, the behavior is dominated by the biggest $|\lambda|$. That is, for large $n$, we must have $(f_{n}, f_{n-1})$ approximately parallel to the corresponding eigenvector, and hence:

$$ \begin{pmatrix} f_{n+1} \\ f_{n} \end{pmatrix} = F \begin{pmatrix} f_{n} \\ f_{n-1} \end{pmatrix} \approx \lambda_1 \begin{pmatrix} f_{n} \\ f_{n-1} \end{pmatrix} $$

where $\lambda_1 = (1 + \sqrt{5})/2$ is the so-called golden ratio.

Let's compute the ratios of $f_{n+1}/f_{n}$ and show that they approach the golden ratio:

In [14]:

(1 + √big(5))/2 # golden ratio computed to many digits

Out[14]:

1.61803398874989484820458683436563811772030917980576286213544862270526046281891

We can also plot the ratio vs. $n$:

In [15]:

fastfib(101) / fastfib(100)

Out[15]:

1.618033988749894848204586834365638117720312743963795685753591851088290198698868

In [16]:

using PyPlot
plot(1:10, [Float64(fastfib(n+1)/fastfib(n)) for n=1:10], "ro-")
plot([0,10], (1+√5)/2 * [1,1], "k--")
xlabel(L"n")
ylabel(L"f_{n+1}/f_n")
title("ratios of successive Fibonacci numbers")

Out[16]:

PyObject Text(0.5, 1.0, 'ratios of successive Fibonacci numbers')

Clearly, it converges rapidly as expected!

(In fact, it converges exponentially rapidly, with the error going exponentially to zero with $n$. We will discuss this in more detail later when discussing the power method.)