I’ve just read yet another article showing how to calculate elements of the Fibonacci sequence in some functional language. I’ll give examples in Scala, though I’m not doing anything very clever linguistically, so they shouldn’t be hard for anyone. Of course, it showed the hopelessly inefficient dual recursion

```
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)
```

and a tail-recursive version

```
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)
```

There was also an explicitly iterative version, which isn’t
interestingly different from the tail-recursive version, so
I shan’t bother showing it. Unfortunately, the tail-recursive
version is still rather slow. On my (very old) laptop, it takes
nearly a minute to calculate `tail_fib(1000000)`

. We can do
*much* better.

## 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.

## Leave a comment