How to really calculate Fibonacci numbers

def terrible_fib(n: BigInt): BigInt =
  if (n < 0) sys.error("not defined for -ve n")
  else if (n == 0) 0
  else if (n == 1) 1
  else terrible_fib(n - 1) + terrible_fib(n - 2)

def tail_fib0(n: BigInt, a: BigInt, b: BigInt): BigInt =
  if (n == 0) a
  else tail_fib0(n - 1, a + b, a)
def tail_fib(n: BigInt): BigInt =
  if (n < 0) sys.error("not defined for -ve n")
  else tail_fib0(n, 0, 1)

Negative elements

Actually, all of the versions in the article I read had a common defect: they all returned incorrect answers for $n \le 0$. So let’s investigate this first.

I’ll define the Fibonacci sequence as $(F_n)_n$, where: \[ F_0 = 0; \qquad F_1 = 1; \qquad F_{n+1} = F_n + F_{n-1} \] In particular, then, $F_0$ is defined and is zero; the first few terms are then $0$, $1$, $1$, $2$, $3$, $5$, $8$, $13$, $21$, $34$, $55$, …

This definition sense for integers $n < 0$: rearranging gives \[ F_{n-1} = F_{n+1} - F_n \] Going backwards from $F_0$, then, gives the terms $0$, $1$, $-1$, $2$, $-3$, $5$, $-8$, … Hey! That looks familiar. Specifically, it’s the same old sequence, only the even items have been negated. Indeed, if we assume inductively that $F_{-n} = (-1)^{n+1} F_n$ for $n \le N$, then \[ \begin{align} F_{-N-1} = F_{1-N} - F_{-N} &= (-1)^N F_{N-1} - (-1)^{N+1} F_N \\ &= (-1)^N (F_{N-1} + F_N) \\ &= (-1)^{N+2} F_{N+1} \end{align} \] as required.

That was the easy part. Now, about actually calculating terms of the Fibonacci sequence efficiently…

Expressing a linear recurrence in matrix form

We start by noticing that each term is a linear function of previous terms. We can express linear functions using matrix multiplication. Therefore, let’s define \[ \mathbf{F} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \] This is interesting because \[ \mathbf{F} \begin{pmatrix} F_n \\ F_{n-1} \end{pmatrix} = \begin{pmatrix} F_n + F_{n-1} \\ F_n \end{pmatrix} = \begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} \] i.e., multiplying a vector containing a consective pair of Fibonacci terms by our matrix $\mathbf{F}$ advances the vector by one step through the sequence. Leaping ahead through the sequence to find some specific future term therefore corresponds to multiplying the vector by an appropriate power of $\mathbf{F}$. And we can calculate powers efficiently using addition chains.

Before we get too excited, there’s more theory to do. We can start by noticing that powers $\mathbf{F}^n$ have a rather special form. Here are the first few. \[ \mathbf{F}^1 = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}; \qquad \mathbf{F}^2 = \begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix}; \qquad \mathbf{F}^3 = \begin{pmatrix} 3 & 2 \\ 2 & 1 \end{pmatrix}; \qquad \mathbf{F}^4 = \begin{pmatrix} 5 & 3 \\ 3 & 2 \end{pmatrix} \] Uhuh… So maybe \[ \mathbf{F} \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix} = \begin{pmatrix} F_{n+1} + F_n & F_n + F_{n-1} \\ F_{n+1} & F_n \end{pmatrix} = \begin{pmatrix} F_{n+2} & F_{n+1} \\ F_{n+1} & F_n \end{pmatrix} \] OK, let’s generalize our matrix a bit, and set \[ \mathbf{F}_n = \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix} \] Then our friend $\mathbf{F} = \mathbf{F}_1$, and, based on what we just calculated, \[ \mathbf{F}_n = \mathbf{F}^n \] We’ve established this for $n \ge 1$ already, and a quick check shows that \[ \mathbf{F}_0 = \begin{pmatrix} F_1 & F_0 \\ F_0 & F_{-1} \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} = \mathbf{I} = \mathbf{F}^0 \] which leaves only one more thing to check: \[ \mathbf{F}_{-1} \mathbf{F}_1 = \begin{pmatrix} F_0 & F_{-1} \\ F_{-1} & F_{-2} \end{pmatrix} \begin{pmatrix} F_2 & F_1 \\ F_1 & F_0 \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & -1 \end{pmatrix} \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} = \mathbf{I} = \mathbf{F}^0 \] and hence $\mathbf{F}_{-1} = \mathbf{F}^{-1}$.

So in fact all we need to do to find the $n$th Fibonacci number is to calculate the top-right element of $\mathbf{F}^n$.

Calculating powers

OK, so how do we calculate powers? There’s a lot of theory here, and it would lead us a long way from the main point, so I’ll just describe one fairly simple approach: right-to-left square-and-multiply.

Suppose we want to calculate $x^n$ for some integer $x$ and some non-negative integer $n$. We start by writing $n$ down as a $w$-bit binary number $(n_{w-1} n_{w-2} \cdots n_1 n_0)_2$, so that $n = \sum_{0\le i< w} n_i 2^i$. Operationally, we start by setting $t_w = 1$. We then loop, for each $i$ from $w - 1$ down to $0$: if $n_i = 0$, then we square, setting $t_i = t_{i+1}^2$; if $n_i = 1$ then we also multiply, setting $t_i = x t_{i+1}^2$. The final result is $t_0 = x^n$.

Maybe you can see intuitively why this works, but it’s always better to have a proof. Pick some bit position $i$ (where $0 \le i < w$). Write $u_i = \sum_{0\le j< i} n_j 2^j$ and $v_i = \sum_{i\le j< w} n_j 2^{j-i}$; then $n = u_i + v_i 2^i$.

We assume, inductively, that $t_{i+1} = x^{v_{i+1}}$. This is certainly true initially, since $v_w = 0$ and $t_w = x^0 = 1$. Observe that $v_i = n_i + 2 v_{i+1}$. Then we set $t_i = x^{n_i} t_{i+1}^2 = x^{n_i} (x^{v_{i+i}})^2 = x^{n_i + 2 v_{i+1}} = x^{v_i}$ as required.

In code, then, we get this.

def pow(x: BigInt, n: BigInt): BigInt = {
  def pow0(i: Int, t: BigInt): BigInt =
    if (i < 0) t
    else if (!n.testBit(i)) pow0(i - 1, t*t)
    else pow0(i - 1, x*t*t);
  pow0(n.bitLength, 1)
}

Which is great if we wanted to calculate powers of integers; but we don’t. We want to calculate powers of a matrix. So we’ll have to factor out the details of integer arithmetic. We only have two — multiplication-by-$x$, and squaring — so we’ll just pass them in as functions. This gives us the following definition.

def genpow[T](n: BigInt, t_init: T, mulx: T => T, sqr: T => T): T = {
  def genpow0(i: Int, t: T): T =
    if (i < 0) t
    else if (!n.testBit(i)) genpow0(i - 1, sqr(t))
    else genpow0(i - 1, mulx(sqr(t)));
  genpow0(n.bitLength, t_init)
}

Putting it all together

The remaining difficult part is the matrix arithmetic. Notice that the matrices $\mathbf{F}^n$ have a special form. We can capture the salient details with just the two numbers in the right-hand column. If we write \[ M(a, b) = \begin{pmatrix} a + b & a \\ a & b \end{pmatrix} \] then the two operations we need, namely multiplication by $\mathbf{F}$ and squaring, work out like this. \[ \begin{gather} \mathbf{F} M(a, b) = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} a + b & a \\ a & b \end{pmatrix} = \begin{pmatrix} 2 a + b & a + b \\ a + b & a \end{pmatrix} = M(a + b, a) \\ M(a, b)^2 = \begin{pmatrix} a + b & a \\ a & b \end{pmatrix}^2 = \begin{pmatrix} 2 a^2 + 2 a b + b^2 & a^2 + 2 a b \\ a^2 + 2 a b & a^2 + b^2 \end{pmatrix} = M(a^2 + 2 a b, a^b + b^2) \end{gather} \] And finally, we end up with the following code.

def genpow[T](n: BigInt, t_init: T, mulx: T => T, sqr: T => T): T = {
  def genpow0(i: Int, t: T): T =
    if (i < 0) t
    else if (!n.testBit(i)) genpow0(i - 1, sqr(t))
    else genpow0(i - 1, mulx(sqr(t)));
  genpow0(n.bitLength, t_init)
}

def fib(n: BigInt): BigInt = {
  val s = if (n < 0 && !n.testBit(0)) -1 else +1;
  val (a, _) = genpow[(BigInt, BigInt)](n.abs, (0, 1),
    { case (a, b) => (a + b, a) },
    { case (a, b) => val aa = a*a; (aa + 2*a*b, aa + b*b) });
  a*s
}

If we don’t care about reuse, we can squish these together into a single function which just does the whole job. This also has the benefit of avoiding messing about with type parameters and lambda- expressions.

def fib(n: BigInt): BigInt = {
  val n1 = n.abs;
  val s = if (n < 0 && !n.testBit(0)) -1 else +1;
  def fib0(i: Int, a: BigInt, b: BigInt): BigInt =
    if (i < 0) a
    else {
      val aa = a*a;
      val a1 = aa + 2*a*b;
      val b1 = aa + b*b;
      if (!n1.testBit(i)) fib0(i - 1, a1, b1)
      else fib0(i - 1, a1 + b1, a1)
    };
  s*fib0(n1.bitLength, 0, 1)
}

And we’re done.

Another approach using polynomials

Oh, before I finish, there’s a slicker way to get to the final formulae, without messing around with matrices. We start by supposing that there’s a number $\varphi$ such that $\varphi^2 = \varphi + 1$. (Such a number exists: if you solve the quadratic, you find that $\varphi = (1 \pm \sqrt{5})/2$, which I hope explains the name.) Let’s consider various powers of $\varphi$ for a moment. \[ \begin{align} \varphi^0 &= 1 = 0 \cdot \varphi + 1; \\ \varphi^1 &= \varphi = 1 \cdot \varphi + 0; \\ \varphi^2 &= \varphi + 1 = 1 \cdot \varphi + 1; \\ \varphi^3 &= \varphi \cdot \varphi^2 = \varphi^2 + \varphi = \varphi + 1 + \varphi = 2 \cdot \varphi + 1; \\ \varphi^4 &= \varphi \cdot \varphi^3 = 2 \varphi^2 + \varphi = 2 \varphi + 2 + \varphi = 3 \cdot \varphi + 2; \\ \varphi^5 &= \varphi \cdot \varphi^4 = 3 \varphi^2 + 2 \varphi = 3 \varphi + 3 + 2 \varphi = 5 \cdot \varphi + 3 \end{align} \] Hmmm. Maybe… \[ \varphi (F_n \varphi + F_{n-1}) = F_n \varphi^2 + F_{n-1} \varphi = F_n (\varphi + 1) + F_{n-1} \varphi = (F_n + F_{n-1}) \varphi + F_n \] And, yes, $\varphi^0 = 1 = F_0 \varphi + F_{-1}$. From the defining equation, $\varphi^2 = \varphi + 1$, so $1 = \varphi^2 - \varphi$ and (since $\varphi \ne 1$), $\varphi^{-1} = \varphi - 1 = F_{-1} \cdot \varphi + F_{-2}$.

We’ve just proven that \[ \varphi^n = F_n \varphi + F_{n-1} \] Hence, if we reduce $\varphi^n$ down to a linear expression in $\varphi$ using the defining equation $\varphi^2 = \varphi + 1$, then the resulting coefficient of $\varphi$ is the $n$th Fibonacci number.

Again, we can go about this by squaring and multiplying, but this time, we’ll work with linear expressions in $\varphi$, i.e., terms of the form $a \varphi + b$. We need to be able to multiply by $\varphi$, and square, so: \[ \begin{gather} \varphi (a \varphi + b) = a \varphi^2 + b \varphi = a (\varphi + 1) + b \varphi = (a + b) \varphi + a \\ (a \varphi + b)^2 = a^2 \varphi^2 + 2 a b \varphi + b^2 = a^2 (\varphi + 1) + 2 a b \varphi + b^2 = (a^2 + 2 a b) \varphi + (a^2 + b^2) \end{gather} \] which leads to exactly the same calculation as we used above.