There was some discussion on reddit recently of the Fibonacci numbers (1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1,597, 2,584, 4,181, 6,765, 10,946, 17,711, 28,657, 46,368, 75,025, 121,393, 196,418, 317,811, 514,229, 832,040, …) and efficient ways of calculating them.

One way of doing so is using numbers of the form *a* + *b* σ where σ is the square root of 5. Multiplication of such numbers satisfies:

(*a* + *b* σ) × (*c* + *d* σ) = *ac* + 5*bd* + (*ad* + *bc*) σ.

We can define the golden ratio φ = (1 + σ) / 2 and also ψ = 1 − φ = (1 − σ) / 2, in which case the *n*^{th} Fibonacci number *F*_{n} will be exactly (φ^{n} − ψ^{n}) / σ. This is known as Binet’s formula.

We can use this formula to calculate Fibonacci numbers using only integer arithmetic, without ever evaluating the σ. We will have:

(2 φ)^{n} − (2 ψ)^{n} = (1 + σ)^{n} − (1 − σ)^{n} = 0 + *p* σ

for some integer *p*, and division by a power of two will give *F*_{n} = *p* / 2^{n}.

I am using the R language, with the `gmp`

package, which provides support for large integer matrices, and this allows us to use the relationship:

If we call this matrix *A* and calculate *A*^{n−1}, the first number in the resultant matrix will be the *n*^{th} Fibonacci number *F*_{n}. The following R code calculates *F*_{10} = 55 using a combination of multiplication and squaring:

```
n <- 10
A <- matrix.bigz(c(
1, 1,
1, 0), 2)
p <- function(n) {
if (n == 1) A
else if (n %% 2 == 1) A %*% p(n-1)
else {
b <- p(n/2)
b %*% b
}
}
p(n-1)[1,1]
```

This same code will calculate, for example:

The time taken to calculate *F*_{n} is approximately proportional to *n*^{1.156}, with the case of *n* = 1,000,000,000 (giving a number with 208,987,640 digits) taking about a minute.