Exhibit: Liber Abacus

There is hardly any sequence of numbers, except maybe the primes which is as famous as the Fibonacci sequence. Similarly, there is hardly any real number besides π and the Euler-e which is as famous as the golden ratio φ = (1+5^1/2)/2. Linear algebra relates the Fibonacci sequence with the golden ratio. If the recursion x(n+1) = x(n) + x(n-1) with the initial condition x(0)=0,x(1)=1 is solved by writing it as a discrete dynamical system v(n+1) = A v(n) written out as

    | x(n+1) |    | 1   1  |   | x(n)   |     | x(1) |   | 1 |
    | x(n)   |  = | 1   0  |   | x(n-1) |,    | x(0) | = | 0 |

The eigenvalues are φ=1.61803... and -1/φ. To solve the system, we write the initial condition as a linear combination of the eigenvectors

| 1 |      |  φ   |         | -1/φ |
| 0 |  = ( |   1  |    -    |   1  |  )/ 5^1/2

then write down the closed form solution

| x(n+1) |        |  φ |              | -1/φ  |
| x(n)   | = (φⁿ  |  1 |  -  (-1/φ)ⁿ  |    1  |  )/ 5^1/2

Looking at the second coordinates gives

 x(n) = [( (1+5^1/2)/2)ⁿ - (1-5^1/2)/2)ⁿ ]/5^1/2

In Mathematica

 phi=(1+Sqrt[5])/2; 
 F[n_]:=(phi^n - (-1/phi)^n)/Sqrt[5];
 Table[Simplify[F[n]],{n,10}]

Now, Mathematica has this already implemented. The same result can be obtained with

  Table[Fibonacci[n],{n,10}]

We could also program the recursion

  G[0]=0; G[1]=1; G[n_]:=G[n-1]+G[n-2];
  Table[G[n],{n,10}]

Now, you might think, whats the point? If the procedure is already built into the system or can be programmed with recursion of half a line, why do we need the explicit solution? Lets look how long a machine needs to compute the built in Fibonacci number for n=1000000000 or the recursion or then using the F formula above:

First[Timing[Fibonacci[10^9]]]

Now, the recursion with G fails with an error "Recursion depth of 1024 exceeded during evaluation". Well, lets increase the recursion limit:

$RecursionLimit=10^9; G[0]=0; G[1]=1; G[n_]:=G[n-1]+G[n-2]; First[Timing[G[10^9]]]

If you run it you see it kills the program. One of the ways to segfault Mathematica! What is nifty about the explicit formula is that we can just say that F[1000000000] is equal to (phi^(1000000000) - (-1/phi)^(100000000))/Sqrt[5]. This takes 50 milli seconds:

phi=(1+Sqrt[5])/2; F[n_]:=(phi^n - (-1/phi)^n)/Sqrt[5]; First[Timing[F[10^9]]]

However, doing the same command with 10^13 gives an error "Uncaught SystemException returned to top level". The memory resources to handle such large algebraic expression was too large.

Folio 139v-140r of the Liber Abbaci. (Biblioteca Nazionale Centrale, Florence, Italy)

Source