In early 2019, due to a silly flamewar, some friends in Taiwan and I took an interest in computation of Fibonacci numbers. This post involving some inductive proofs and some light program derivation. If you think the fastest way to compute Fibonacci numbers is by a closed-form formula, you should read on!

Let `Nat`

be the type of natural numbers. We shall all be familiar with the following definition of Fibonacci number:

```
``` fib :: Nat -> Nat
fib 0 = 0
fib 1 = 1
fib (n+2) = fib (n+1) + fib n

`(When defining functions on natural numbers I prefer to see `

`0`

, `(+1)`

(and thus `(+2)`

), as constructors that can appear on the LHS, while avoiding subtraction on the RHS. It makes some proofs more natural, and it is not hard to recover the Haskell definition anyway.)

= (+1) . (+1)

Executing the definition without other support (such as memoization) gives you a very slow algorithm, due to lots of re-computation. I had some programming textbooks in the 80’s wrongly using this as an evidence that “recursion is slow” (`fib`

is usually one of the only two examples in a sole chapter on recursion in such books, the other being tree traversal).

By defining `fib2 n = (fib n, fib (n+1))`

, one can easily derive an inductive definition of `fib2`

,

```
``` fib2 :: Nat -> (Nat, Nat)
fib2 0 = (0, 1)
fib2 (n+1) = (y, x+y)
where (x,y) = fib2 n

` which computes `

`fib n`

(and `fib (n+1)`

) in `O(n)`

recursive calls. Be warned, however, that it does not imply that `fib2 n`

runs in `O(n)`

time, as we shall see soon.

## Binet’s Formula

To be even faster, some might recall, do we not have a closed-form formula for Fibonacci numbers?

```
``` fib n = (((1+√5)/2)^n - ((1-√5)/2)^n) /√5

`It was believed that the formula was discovered by Jacques P. M. Binet in 1843, thus we call it `

*Binet’s formula* by convention, although the formula can be traced back earlier. Proving (or even discovering) the formula is a very good exercise in inductive proofs. On that I recommend this tutorial by Joe Halpern (CS 280 @ Cornell, , 2005). Having a closed-form formula gives one an impression that it give you a quick algorithm. Some even claim that it delivers a `O(1)`

algorithm for computing Fibonacci numbers. One shall not assume, however, that `((1+√5)/2)^n`

and `((1-√5)/2)^n`

can always be computed in a snap!

When processing large numbers, we cannot assume that arithmetic operations such as addition and multiplication take constant time. In fact, it is fascinating knowing that multiplying large numbers, something that appears to be the most fundamental, is a research topic that can still see new breakthrough in 2019 [HvdH19].

## Vorobev’s Equation

There is another family of algorithms that manages to compute `fib n`

in `O(log n)`

recursive calls. To construct such algorithms, one might start by asking oneself: can we express `fib (n+k)`

in terms of `fib n`

and `fib k`

(and some other nearby `fib`

if necessary)? Given such a formula, we can perhaps compute `fib (n+n)`

from `fib n`

, and design an algorithm that uses only `O(log n)`

recursive calls.

Indeed, for `n >= 1`

, we have

```
``` fib (n+k) = fib (n-1) * fib k + fib n * fib (k+1) . -- (Vor)

`This property can be traced back to Nikolai. N. Vorobev, and we therefore refer to it as `

*Vorobev’s Equation*. A proof will be given later. For now, let us see how it helps us.

With Vorobev’s equation we can derive a number of (similar) algorithms that computes `fib n`

in `O(log n)`

recursive calls. For example, let `n, k`

in (Vor) be `n+1, n`

, we get

```
``` fib (2n+1) = (fib (n+1))^2 + (fib n)^2 -- (1)

`Let `

`n, k`

be `n+1, n+1`

, we get

```
``` fib (2n+2) = 2 * fib n * fib (n+1) + (fib (n+1))^2 -- (2)

`Subtract (1) from (2), we get`

```
``` fib 2n = 2 * fib n * fib (n+1) - (fib n)^2 -- (3)

The LHS of (1) and (3) are respectively odd and even, while their RHSs involve only `fib n`

and `fib (n+1)`

. Define `fib2v n = (fib n, fib (n+1))`

, we can derive the program below, which uses only `O(log n)`

recursive calls.

```
``` fib2v :: Nat -> (Nat, Nat)
fib2v 0 = (0, 1)
fib2v n | n `mod` 2 == 0 = (c,d)
| otherwise = (d, c + d)
where (a, b) = fib2v (div n 2)
c = 2 * a * b - a * a
d = a * a + b * b

## Which Runs Faster?

Having so many algorithms, the ultimate question is: which runs faster?

Interestingly, in 1988, James L. Holloway devoted an entire Master’s thesis to analysis and benchmarking of algorithms computing Fibonacci numbers. The thesis reviewed algorithms including (counterparts of) all those mentioned in this post so far, and some more algorithms based on matrix multiplication. I will summarise some of his results below.

For a theoretical analysis, we need know the number of bits needed to represent `fib n`

. Holloway estimated that to represent `fib n`

, we need approximately `n * 0.69424`

bits. We will denote this number by `N n`

. That `N n`

is linear in `n`

is consistent with our impression that `fib n`

grows exponentially in `n`

.

Algorithm `fib2`

makes `O(n)`

recursive calls, but it does not mean that the running time is `O(n)`

. Instead, `fib2 n`

needs around `N (n^2/2 - n/2)`

bit operations to compute. (Note that we are not talking about big-O here, but an approximated upper bound.)

What about Binet formula? We can compute `√5`

by Newton’s method. One can assume that each `n`

bit division needs `n^2`

operations. In each round, however, we need only the most significant `N n + log n`

bits. Overall, the number of bit operations needed to compute Binet formula is dominated by `log n * (N n + log n)^2`

— not faster than `fib2`

.

Holloway studied several matrix based algorithm. Generally, they need around `(N n)^2`

bit operations, multiplied by different constants.

Meanwhile, algorithms based on Vorobev’s Equation perform quite well — it takes about `1/2 * (N n)^2`

bit operations to compute `fib2v n`

!

What about benchmarking? Holloway ran each algorithm up to five minutes. In one of the experiments, the program based on Binet’s formula exceeds 5 minutes when `log n = 7`

. The program based on `fib2`

terminated within 5 minutes until `log n = 15`

. In another experiment (using simpler programs considering only cases when `n`

is a power of `2`

), the program based on Binet’s formula exceeds 5 minutes when `log n = 13`

. Meanwhile the matrix-based algorithms terminated within 3 to 5 seconds, while the program based on Vorobev’s Equation terminated within around 2 seconds.

## Proving Vorobev’s Equation

Finally, let us see how Vorobev’s Equation can be proved. We perform induction on `n`

. The cases when `n := 1`

and `2`

can be easily established. Assuming the equation holds for `n`

(that is, (Vor)) and `n:= n+1`

(abbreviating `fib`

to `f`

):

```
``` f (n+1+k) = f n * f k + f (n+1) * f(k+1) -- (Vor')

`we prove the case for `

`n:=n+2`

:

```
``` f (n+2+k)
= { definition of f }
f (n+k) + f (n+k+1)
= { (Vor) & (Vor') }
f (n-1) * f k + f n * f (k+1) +
f n * f k + f (n+1) * f(k+1)
= { f (n+1) = f n + f (n-1) }
f (n+1) * f k + f n * f (k+1) + f (n+1) * f (k+1)
= { f (n+2) = f (n+1) + f n }
f (n+1) * f k + f (n+2) * f (k+1) .

`Thus completes the proof.`

## Related Work

Dijkstra derived another algorithm that computes `fib n`

in `O(log n)`

recursive calls in EWD654 [Dij78].

Besides his master’s thesis, Holloway and his supervisor Paul Cull also published a journal version of their results [CH89]. I do not know the whereabouts of Holloway — it seems that he didn’t pursue a career in academics. I wish him all the best. It comforts me imagining that any thesis that is written with enthusiasm and love, whatever the topic, will eventually be found by some readers who are also enthusiastic about it, somewhere, sometime.

I found many interesting information on this page hosted by Ron Knott from University of Surrey, and would recommend it too.

After the flamewar, Yoda Lee (李祐棠) conducted many experiments computing Fibonacci, taking into considerations things like precision of floating point computation and choosing suitable floating point libraries. It is worth a read too. (In Chinese.)

So, what was the flamewar about? It started with someone suggesting that we should store on the moon (yes, the moon. Don’t ask me why) some important constants such as `π`

and `e`

and, with the constants being available in very large precision, many problems can be solved in constant time. Then people started arguing what it means computing something in constant time, whether Binet’s formula gives you a constant time algorithm… and here we are. Silly, but we learned something fun.

## References

[**CH89**] Paul Cull, James L. Holloway. Computing fibonacci numbers quickly. Information Processing Letters, 32(3), pp 143-149. 1989.

[**Dij78**] Dijkstra. In honor of Fibonacci. EWD654, 1978.

[**Hol88**] James L. Holloway. Algorithms for Computing Fibonaci Numbers Quickly. Master Thesis, Oregon State University, 1988.

[**HvdH19**] David Harvey, Joris Van Der Hoeven. Integer multiplication in time `O(n log n)`

. 2019. hal-02070778.