A novel and efficient way to compute Fibonacci numbers

An earlier post described how to compute Fibonacci numbers in a single arithmetic expression.

Faré Rideau, the author of a page of Fibonacci computations in Lisp, suggested in a private email a simple and efficient variant, that I believe is novel.

For $X$ large enough, $\mathrm{Fib}_n = (X^{n+1}\ \mathrm{mod}\ (X^2-X-1))\ \mathrm{mod}\ X$.

That means you can compute Fibonacci numbers efficiently with a simple program:

for n in range(1, 21):
    X = 1<<(n+2)
    print(pow(X, n+1, X*X-X-1) % X)

This blog post describes how this method works, gives a few ways to think about it, easily infers the fast Fibonacci doubling formulas, provides a nice alternative to Binet’s formula relating the golden ratio and Fibonacci numbers, and expands the method to generalized Fibonacci recurrences, including a near one-line solution to the problem of counting how many ways to reach the end-square of a 100-square game using a single six-sided dice.

Computing with integers

First, if you’ve not read my earlier Fibonacci article then go ahead and read it now, especially the overview which covers the common ways to generate Fibonacci numbers.

To see how Faré’s formula works, let’s see what happens if we calculate $(aX + b)X$ modulo $X^2-X-1$. Modulo $X^2-X-1$, $X^2 = X+1$, so we can replace any $X^2$ that appear by $X+1$: $$ (aX + b)X = aX^2 + bX = a(X+1) + bX = (a+b)X + a $$

If $a$ and $b$ are small relative to $X$ – specifically if $(a+b)X + b < X^2 - X - 1$ – then we’ve reproduced one step of the standard iterative method for computing Fibonacci: $a, b \rightarrow a+b, a$. (If $a$ and $b$ get too large relative to $X$, then the equations still hold, but some of the $(a+b)X$ term will overflow, wrapping around into the $a$ term).

If we set $a=\mathrm{Fib}_1=1, b=\mathrm{Fib}_0=0$, so that $aX+b = X$ we can apply our calculation above $n$ times to get $X^n = \mathrm{Fib}_nX + \mathrm{Fib}_{n-1}\ \mathrm{mod}\ X^2-X-1$.

As a worked example of this method, let’s see how it looks with $X=100$. Then $X^2-X-1 = 9899$. Tabulating powers of $100$ and reducing modulo $9899$:

We can see in the final column that $\mathrm{Fib}_{n-1}$ appears in the lower two digits (0, 1, 1, 2, 3, 5, …, 34), and $\mathrm{Fib}_n$ appears in the upper two digits (1, 1, 2, 3, 5, 8, …, 55).

So, If we want to compute $\mathrm{Fib}_n$, we need to pick $X$ large enough (so that $\mathrm{Fib}_nX +\ \mathrm{Fib}_{n-1} < X^2-X-1$, compute $X^{n+1}\ \mathrm{mod}\ X^2-X-1$, and take the result modulo $X$. Faré uses a tight bound for $X$, but $X=2^{n+2}$ is easier (at the cost of a little efficiency), and that’s what I’ll use here.

In Python:

def fib(n):
	X = 1<<(n+2)
    return pow(X, n+1, X*X-X-1) % X

for i in range(1, 21):
    print(i, fib(i))        

This code uses Python’s built-in ternary pow(a, b, m) function that efficiently computes the result of a to the power of b modulo m in exact integer arithmetic, using exponentiation by squaring.

Computation here is done using approximately $2n$ bits, which is reasonable since Fibonacci numbers grow exponentially, so the number of bits grows linearly.

I find this code quite beautiful. Not only is it efficient and relatively simple to read (if not to understand), it allows one to compute Fibonacci numbers with no (explicit) iteration or recursion.

Computing in polynomials

We can also compute the result using this method symbolically, treating $X$ as an unknown. Then our “numbers” are linear polynomials of the form $a + bX$, and multiplication is polynomial multiplication modulo $X^2-X-1$. If you’re more mathematically sophisticated, this corresponds to the quotient ring $\mathbb{Z}[X] / (X^2-X-1)$.

We can calculate what multiplication looks like here: $(a_0 + a_1X)(b_0+b_1X) = a_0b_0 + (a_1b_0+a_0b_1)X + a_1b_1X^2$. As before, using that $X^2=X+1$, this equals $(a_0b_0+a_1b_1) + (a_1b_0+a_0b_1+a_1b_1)X$.

We can write some code that implements this multiplication, and performs exponentiation by squaring using it.

def mul(a, b):
	return a[0]*b[0]+a[1]*b[1], a[1]*b[0]+a[0]*b[1]+a[1]*b[1]

def ppow(a, n):
	r = (1, 0)
	while n:
        if n & 1: r = mul(r, a)
        a = mul(a, a)
        n >>= 1
	return r

for i in range(1, 21):
	print(i, ppow((0, 1), i+1)[0])

This code is very like the matrix-power method for computing Fibonacci numbers, and it’s interesting to think of how the two correspond.

We can also reproduce the fast doubling formulas, which are another good method for computing Fibonacci numbers. The fast doubling formulas are: $$ \begin{align*} & \mathrm{Fib}_{2n} = \mathrm{Fib}_n(2\mathrm{Fib}_{n+1} - \mathrm{Fib}_n) \\ & \mathrm{Fib}_{2n+1} = \mathrm{Fib}_{n+1}^2 + \mathrm{Fib}_n^2\end{align*}$$

These drop out straightforwardly from calculation using our polynomial method. With all calculations being performed modulo $X^2-X-1$ we have: $\mathrm{Fib}_{2n}X + \mathrm{Fib}_{2n-1} = X^{2n} = X^nX^n = (\mathrm{Fib}_nX + \mathrm{Fib}_{n-1})^2$. Expanding out this last term and reducing it modulo $X^2-X-1$ gives us something that’s easily massaged into the doubling formulas.

Computing with the golden ratio

Another (less practical) way to think of this method is to consider numbers of the form $a + b\phi$, where $\phi$ is the golden ratio $\frac{1+\sqrt{5}}{2}$. Since $\phi^2 = \phi + 1$, multiplications in this system works out identically to how it did in the linearly polynomial case, simply replacing $X$ with $\phi$. For the more mathematically sophisticated, that’s noting that $\mathbb{Q}[X] / (X^2-X-1) \simeq \mathbb{Q}[\phi]$.

A nice alternative to Binet’s formula that directly relates the golden ratio with Fibonacci numbers that comes from this is $\phi^n = \mathrm{Fib}_n\phi + \mathrm{Fib}_{n-1}$.

Generalizing to $n$-acci numbers

These methods can be straightforwardly generalized to higher-order Fibonacci numbers.

For example, suppose we want to calculate the number of ways to get from the start to (exactly) square 100 of a snakes-and-ladders game (with no snakes or ladders) using a single dice, we have the recurrences:

R(n) = 0 (n<0)
R(0) = 1
R(n) = R(n-1) + R(n-2) + R(n-3) + R(n-4) + R(n-5) + R(n-6) for n > 1

These are the hideously named Hexanacci numbers (although they’re usually defined as a sequence starting 0, 0, 0, 0, 0, 1, which is offset from 5 from R). See the Wikipedia article on higher order Fibonacci numbers and sequence AA01592 on OEIS.

We can use the same tricks as before to calculate this sequence in very little code, and quite efficiently, using $log(n)$ multiplications of numbers with approximately $6n$ bits. Again, we can’t expect much better than this, since the results grow exponentially.

def hexa(n):
	X = 1<<(n+1)
	return pow(X, n+6, X**6-X**5-X**4-X**3-X**2-X-1) % X

for i in range(1, 21):
	print(i, hexa(i))

Computing hexa(100) with this function gives 290078479914610587823630044098.

If we wanted, we could easily turn this computation into a single line program.

Further thinking required

Here’s some (as far as I know open) and interesting questions.

• How general is this approach? What other recurrences can be computed like this?
• How does one compute Fibonacci-like numbers with different starting values? Can you write down a formula for the $n$-acci numbers starting with an arbitrary $n$ values $a_0, a_1, …, a_n$?
• The fast doubling formulas for the Fibonacci numbers drop out of the polynomial power method. Can you write down fast doubling formulas for $n$-acci numbers?
• How do you efficiently compute the polynomial multiplications for the $n$-acci numbers?
• Compare the cost (in the bit-cost model) of the matrix power method, integer power method, and polynomial power method for computing $n$-acci numbers. Which is the most efficient?