Calculating The Billionth Fibonacci Number

After reading the title, many of you might wonder why anyone would take on such a task. Well, because I can. Sometimes, you just have to challenge yourself to do something ambitious. While it may sound straightforward, there are actually many challenges to consider. You can't just use any arbitrary algorithm and let it run for a billion iterations, depending on the method, that could take years, maybe even decades. Another issue is storing the number. The billionth Fibonacci number is so incredibly large that it's beyond human comprehension, and even computers struggle to handle it. So, I wanted to see if I could calculate it. But not only that, I also wanted it to be fast. So, let's dive right in.

First of all, let's talk a bit more about how the Fibonacci numbers work. For those not already familiar with the concept, it is quite simple. Each Fibonacci number is calculated by summing the previous two Fibonacci numbers. From now on, we will define $\text{Fib}(n)$ to mean the $n^{\text{th}}$ Fibonacci number. Mathematically, we can write this as follows:

\text{Fib}(n) = \text{Fib}(n-1) + \text{Fib}(n-2)

The following equations are the base cases for the Fibonacci sequence. These are always the same:

\text{Fib}(0) = 0 \quad \text{and} \quad \text{Fib}(1) = 1

Now, the calculations are quite simple. However, as we enter the billion range, even such a simple equation can become prohibitively expensive to calculate. There's also the little problem of computer data types. Most of the times they are limited to a certain size. For example, your standard integer is just 32-bit long. That's barely enough to hold the number 4 billion. Since Fibonacci numbers grow in size at an exponential rate, our traditional data types unfortunately won't work for this endeavor.

There are multiple ways to tackle this problem. The first would be coding your own arbitrary-sized data types. However, that's quite a headache, and I'm lazy. We could also use a little library called GMP for this (more on this later); however, for now, we're going to go even simpler. The perfect solution here is Python. Yes, our green little snakey friend. Python natively supports arbitrary-sized integers. So, just using Python's standard int type will work for our program. I will also take this as an opportunity to prove that while yes, Python might be slower than many languages, with the right modifications, it can be more than adequately fast for most endeavors.

Before we start, the goal is not only to calculate the Billionth Fibonacci Number, but rather to calculate it in a relatively short timespan. What exactly that time span is will be talked about more later. However I can tell you that I dont wanna have to wait years for the result. Another important detail: All the benchmarks and performance data you will see here were performed on my M2 Pro MacBook Pro with 16GB of RAM. With that out of the way let's now look at a few approaches to calculate the Fibonacci sequence.

Recursive

The first and simplest solution is simply writing the code exactly as the equation. However, this is also by far the worst solution out there. But, let's first look at the code:

pythondef fib(n):
    if n <= 1:
        return n
    return fib(n-1) + fib(n-2)

It is a very simple and beautiful recursive function. However, as it stands right now, in many cases it isn't even a solution at all. For example, for $n = 43$ , the program already takes over 40 seconds to complete. But why is that? Well, for each function call, it calls itself twice, where each call then again calls itself twice. This leads to exponential growth where, in total, we have to perform $2^n$ operations. In mathematical terms, it has a runtime complexity of $O(2^n)$ . Meaning, for an input of 43, the method called itself $8.796.093.022.208$ times. You can imagine what $2^{1.000.000.000}$ would look like. Since 43 is a bit shy of one billion, this method won't work for us.

Performance of Recursive Method

Iterative

A bit more efficient but still relatively easy to understand. The iterative approach eliminates one of the most critical bottlenecks of the recursive approach by defining $\text{fib}(0)$ and $\text{fib}(1)$ at the start and then simply keeping track of the last two calculated Fibonacci numbers and summing them up. Here's the code for this solution:

pythondef fib(n):
    a, b = 0, 1  # Initialize the first two Fibonacci numbers
    for _ in range(n):
        a, b = b, a + b
    return a

And yes, this time it is an actual solution. While still being incredibly slow, it would reach the billionth Fibonacci number eventually. Big emphasis on "Eventually." Compared to the recursive solution, however, it sports an incredibly runtime complexity of "only" $O(n)$ . While not great, it's still a massive improvement over the recursive method and would only need a billion total iterations to reach our goal. But since I don't quite have time to wait for a billion iterations either, I will still consider this method a fail.

Performance of Iterative Method

Matrix Exponentiation

Now this is where things start to get interesting. Instead of looking at this problem on a high level, we are now starting to look a bit lower. One property of the Fibonacci sequence that will be very beneficial to us is the fact that it can be calculated via matrix exponentiation. We can represent this as:

\begin{bmatrix}1 & 1 \\1 & 0\end{bmatrix}^n = \begin{bmatrix}F_{n+1} & F_n \\F_n & F_{n-1}\end{bmatrix}

Based on this, we can calculate the $n^{\text{th}}$ Fibonacci number by taking the matrix to the $n^{\text{th}}$ power. The good thing here is that computers are fast at matrix operations. Now, naively, matrix exponentiation will still take up $O(n)$ matrix multiplications. Now if we factor in that each single matrix multiplication does multiple additions and multiplications we actually are going to end up with more total calculations than the iterative approach. However, with a few clever tricks, you can calculate matrix exponentiation in $O(\log n)$ time.

pythondef fibonacci_matrix(n):
    F = [[1, 1], [1, 0]]  # Base matrix for Fibonacci sequence
    result = matrix_pow(F, n)
    return result[0][1]  # The nth Fibonacci number is in this position
@staticmethod
def matrix_mult(A, B):
    return [
        [A[0][0] * B[0][0] + A[0][1] * B[1][0], A[0][0] * B[0][1] +
            A[0][1] * B[1][1]],
        [A[1][0] * B[0][0] + A[1][1] * B[1][0], A[1][0] * B[0][1] +
            A[1][1] * B[1][1]]
    ]
def matrix_pow(self, mat, p):
    result = [[1, 0], [0, 1]]  # Identity matrix
    while p > 0:
        if p & 1:  # If p is odd
            result = matrix_mult(result, mat)
        mat = matrix_mult(mat, mat)  # Square the matrix
        p >>= 1  # Divide p by 2
    return result

This is called Fast Exponentiation. By implementing this as our exponentiation algorithm, we can get the $n^{\text{th}}$ Fibonacci number in $O(\log n)$ runtime. This is significantly faster than all previous attempts, and it would actually get us our billionth number in a somewhat reasonable time.

Performance of Matrix Exponentiation Method

However, we can do even better...

Improved Matrix Exponentiation

This method, similarly to above, is still based on matrix exponentiation. However, we will use a few key improvements to increase the speed of those exponentiations even further. First of all, the matrix is defined a bit differently this time. While the form might be different, it mathematically still works the same.

\begin{bmatrix}0 & 1 \\1 & 1\end{bmatrix}^n = \begin{bmatrix}F_{n-1} & F_n \\F_n & F_{n+1}\end{bmatrix}

pythondef mat_fib(n):
    a, b, c, d = 0, 1, 1, 1  # Initial matrix [[0, 1], [1, 1]]
    x, z = 0, 1  # Initial vector [0, 1]
    while n:
        n, r = divmod(n, 2)
        if r:  # If n is odd
            x, z = a*x + b*z, c*x + d*z  # Multiply current matrix with vector
        if n:  # If n is not zero
            # Square the current matrix
            a, b, c, d = a*a + b*c, b*(a + d), c*(a + d), c*b + d*d
    return x  # x holds the nth Fibonacci number

The key difference for this method is performing operations on single elements instead of the entire matrix. This, simply put, reduces computational overhead and memory usage, and should increase cache efficiency. While the runtime complexity is still $O(\log n)$ , the total runtime decreases as our program does the same number of iterations but with fewer calculations per iteration, and better cache performance, i.e., becoming faster.

Performance Comparison of Normal and Improved Matrix Methods

Doubling Fibonacci

But what if I told you that we can reduce the number of calculations even further? Of course, this method also starts with some mathematical definitions. First of all, we will define two identities from which we can calculate our sequence.

Double Index Identity for Even Numbers:

F(2k) = F(k) \times \left[ 2F(k+1) - F(k) \right]

Double Index Identity for Odd Numbers:

F(2k+1) = F(k+1)^2 + F(k)^2

This basically lets us break down a problem into two smaller subproblems. Essentially, we are halving the index and thereby halving the required iterations on each iteration. Say we want to calculate $\text{Fib}(100)$ , we could write that as:

F(100) = F(50) \times \left[ 2F(51) - F(50) \right]

Meaning we don't have to calculate any numbers between 100 and 51. Now, this method still requires $O(\log n)$ iterations in total. However, it requires fewer operations per iteration compared to both previous methods, as it ditches the matrix operations completely. As a result, this method is overall more efficient.

The code for this looks like this:

pythondef fibonacci_doubling(n):
    def _fibonacci_doubling(n):
        if n == 0:
            return (0, 1)
        else:
            a, b = _fibonacci_doubling(n >> 1)
            c = a * (2 * b - a)
            d = a * a + b * b
            if n & 1:
                return (d, c + d)
            else:
                return (c, d)
    return _fibonacci_doubling(n)[0]

A little side note here. As some of you might have noticed this method works recursively. However contrary to our first approach it is incredibly efficient as it only calculates each number once and skips many more.

Performance Comparison of Doubleing and Improved Matrix Methods

And this is it. This is the method we will be using to get to the billionth Fibonacci number. It is by far the fastest method we have. However it is still not quite fast enough for my liking. And there's one ace still left up my sleeve. While up until this point we have used Python's native arbitrary-sized integers, these are, as you might have guessed, slow. To mitigate this, we will use a cool library called GMP. GMP allows us to work with arbitrary-sized integers as well. However, GMP was originally written in pure C and has a huge number of mathematical and algorithmic optimizations. Rewriting the code to use this library is rather simple. It just requires a few more declarations in each class, but otherwise, the code remains largely unchanged. Now, let's see if it's any faster.

Performance Comparison of GMP and non-GMP Doubling Methods

And wow. GMP completely blows Python out of the water. For very low $n$ 's, it is slower due to the extra overhead of creating these MPZ objects, but it makes up for this significantly on larger $n$ 's. Now, we can calculate the billionth Fibonacci number in just ~7.38 seconds. What will I do with this? Not quite sure yet. However, it was cool seeing how a relatively simple problem like Fibonacci can be solved in so many different ways, ranging from simple to complex mathematics. It was also interesting seeing how fast Python can be if you work with it right. Yes, a solution written entirely in C would probably be even faster, but it would've also taken me 10x more time to code it. Also, just for fun, I ported each approach to GMP. Below you can find the graph with all algorithms compared in log scale.

Performance Comparison of All GMP Methods (log)

Further Ideas

Now, while I do consider this project a success, I have a few things in mind that could be interesting to try out for an even faster program:

Multithreading for Doubling Fibonacci: You could potentially use multithreading since you have basically two problems running at the same time. Will the overhead be worth it? Will it even speed up at all? All questions for a possible revisit...
Leveraging GPUs: Could this potentially be sped up by using GPUs? What about utilizing tensor cores, which are specifically designed for matrix multiplications? Would it be faster than on the CPU? Could you further parallelize it?
Further Mathematical Optimizations: Are there even faster ways to calculate the Fibonacci sequence? Could you use different mathematical identities to further reduce the number of calculations?

GitHub

In case you wanna further explore this project, implement your own methods to compare, or anything else, feel free to check out all the code shown here, and more, on my GitHub Page

Skills Improved

Python +35

C +25