Generating the Fibonacci sequence using matrices

Tue, 25 Mar 2025

Concepts

Before we start, we need to remember some concepts. First of all, what is a sequence? In mathematics, a sequence is a function whose domain is the set of natural numbers. The elements of the sequence are called terms. For example, the sequence of natural numbers is defined by the function $f (n) = n$ , where $n$ is a natural number. The first term of this sequence is $f (1) = 1$ , the second term is $f (2) = 2$ , and so on. This is not the same as a set, because the order of the elements matters in a sequence, while in a set, the order does not matter. For example, the set of natural numbers is $N = {1, 2, 3, 4, \dots}$ , the order of the elements does not matter, we are only interested in the elements of the set, which contains all the natural numbers. On the other hand, the sequence of natural numbers is $1, 2, 3, 4, \dots, a_{n}$ , the order of the elements matters, and the sequence is always infinite.

The Fibonacci sequence is a recursive sequence defined by the following function:

F_{n} = {1, F_{n - 1} + F_{n - 2}, if n = 1 or n = 2 if n > 2

Sometimes, you may see the Fibonacci sequence defined as:

F_{0} = 1, F_{1} = 1, F_{n} = F_{n - 1} + F_{n - 2}

They are all the same, but the second definition is the most common one. If we substitute the values of $n$ in the function, we will have this sequence:

1, 1, 2, 3, 5, 8, 13, 21, 34, \dots, F_{n}

Generating the Fibonacci sequence

We'll see some ways to generate the Fibonacci sequence, from the worst to the best.

Using recursion

The most common way (and worst in terms of performance) to generate the Fibonacci sequence is using recursion. The recursive function is defined by the same function we saw in the equation above, but implemented in a programming. Here is a C++ implementation of this function:

int fib(int n) {
    if (n <= 2)
        return 1;

    return fib(n - 1) + fib(n - 2);
}

Don't ever use this implementation. It's literally the worst way to generate the Fibonacci sequence. The time complexity of this function is $O (2^{n})$ , which is exponential. This means that the time it takes to calculate the $n$ -th term of the sequence grows exponentially with $n$ . This is because the function calls itself multiple times for the same values, which makes the function calculate the same values multiple times. Very inefficient.

Using dynamic programming

A better way to generate the Fibonacci sequence is using dynamic programming. The idea is to store the values of the terms that have already been calculated in an array and use these values to calculate the next terms, we call this memoization. Here is a C++ implementation of this function:

int fib(int n) {
    int memo[n + 1];

    memo[1] = 1;
    memo[2] = 1;

    for (int i = 3; i <= n; i++)
        memo[i] = memo[i - 1] + memo[i - 2];

    return memo[n];
}

This implementation has a time complexity of $O (n)$ , which is linear. This is because the function calculates each term of the sequence only once and stores the values in an array, so the function doesn't need to calculate the same values multiple times. Much more efficient than the recursive implementation, but the problem with this approach is that the space complexity is $O (n)$ . The array that stores the values always grows with $n$ , which can be a problem if $n$ is very large.

Using iterations

One of the best ways to generate the Fibonacci sequence is using iterations. The idea is to calculate the terms of the sequence one by one, without storing the values in an array. Here is a C++ implementation of this function:

#include <algorithm>

using namespace std;

int fib(int n) {
    int a = 1, b = 1;

    for (int i = 3; i <= n; i++) {
        swap(a, b);
        b += a;
    }

    return b;
}

This implementation has a time complexity of $O (n)$ , which is linear, just like the dynamic programming implementation, but the space complexity is $O (1)$ , which is constant. This is because the function doesn't store the values of the terms in an array, it only stores the last two terms, which makes the function very efficient in terms of space complexity.

Theory

Now that we know how to generate the Fibonacci sequence, let's see a different way to do it. We can use matrices to generate the Fibonacci sequence. The idea is to represent the sequence as a matrix and use matrix multiplication to calculate the terms of the sequence. Note that in the last implementation we used iterations to generate the sequence, where we stored only the last two terms. The new values of the terms were calculated based on the last two terms, which means that the terms of the sequence are dependent on the last two terms, i.e. the terms are linearly dependent:

int new_curr = curr + prev;
int new_prev = curr;

// We can express this same idea using linearly dependent equations:

int new_curr = 1 * curr + 1 * prev;
int new_prev = 1 * curr + 0 * prev;

We can rewrite these equations as a matrix multiplication:

[new_curr new_prev] = [1110] [curr prev]

The important part that interests us is that we can multiply the matrix by itself to calculate the next terms of the sequence, let's substitute the values of the matrix with 13 and 21 (the prev and curr values of the sequence, respectively):

[new_curr new_prev] = [1110] [2113] = [1 \cdot 21 + 1 \cdot 13 1 \cdot 21 + 0 \cdot 13] = [3421]

This means that curr is now 34 and prev is now 21. If we multiply the matrix by itself again, we will get the next terms of the sequence:

[new_curr new_prev] = [1110] [3421] = [1 \cdot 34 + 1 \cdot 21 1 \cdot 34 + 0 \cdot 21] = [5534]

Now curr is 55 and prev is 34, and so on. As matrices are associative, we can multiply the matrix by itself $n$ times to get the $n$ -th term of the sequence from the prev and curr values.

[new_curr new_prev] = [1110]^{n} [curr prev]

If we substitute $n$ by some values, you'll see that the linear multiplication of the matrix by itself will also generate the terms of the Fibonacci sequence. Trying with $n = 2$ :

[1110]^{2} [curr prev] = ([1110] [1110]) [curr prev] = ([1 \cdot 1 + 1 \cdot 1 1 \cdot 1 + 0 \cdot 1 1 \cdot 1 + 1 \cdot 0 1 \cdot 1 + 0 \cdot 0]) [curr prev] = [2111] [curr prev]

If we try with $n = 5$ :

[1110]^{5} [curr prev] = ([1110]^{2} [1110]^{2} [1110]) [curr prev] = ([2111] [2111] [1110]) [curr prev] = [5332] [curr prev]

As you can see, the linear dependent matrix itself gives us the next terms of the sequence. What if we get rid of the curr and prev values? We can do that too. The idea is to multiply the matrix by itself $n$ times, and the result still gives us the next terms of the sequence. The only difference is that we will have to multiply the matrix by itself instead of multiplying the matrix by the curr and prev values.

Instead of multiplying the matrix by itself $n$ times (which the complexity would be $O (n)$ ), we can use the binary exponentiation. The time complexity of this approach is $O (lo g n)$ , which is logarithmic, much better than the previous implementations. As we use just one 2x2 matrix, the space complexity is also constant, $O (1)$ , just like the iterations implementation. This can be represented by the following formula:

[1110]^{n} = [F_{n + 1} F_{n} F_{n} F_{n - 1}]

The equivalent code in C++ would be:

#include <algorithm>

using namespace std;

template <typename T>
struct matrix {
    T m[2][2] = {{1, 1}, {1, 0}};

    matrix friend operator*(const matrix& a, const matrix& b) {
        matrix res;

        for (int i = 0; i < 2; i++) {
            for (int j = 0; j < 2; j++) {
                res.m[i][j] = 0;

                for (int k = 0; k < 2; k++)
                    res.m[i][j] += a.m[i][k] * b.m[k][j];
            }
        }

        return res;
    }

    matrix power(T n) {
        if (n == 1)
            return *this;

        matrix ans { .m = {{1, 0}, {0, 1}} };

        while (n) {
            if (n & 1)
                ans = ans * (*this);

            *this = *this * (*this);
            n >>= 1;
        }

        return ans;
    }
};

template <typename T>
T fib(T n) {
    if (n <= 2)
        return 1;

    matrix<T> res;
    res = res.power(n);

    return res.m[0][1];
}

Using arbitrary precision integers

One problem with this implementation is that we can only get the 93rd term of the Fibonacci sequence using a 64-bit integer (unsigned long long), because the 94th term is greater than $2^{63}$ bits, so it'll overflow. To fix this, we need an arbitrary precision integer. The C++ standard library does not have an arbitrary precision integer, but we have a GNU library called GMP that implements this. We get rid of the templates and use the mpz_class class from the library. The code below is a C++ implementation of the Fibonacci sequence using the GMP library:

#include <gmp.h>
#include <gmpxx.h>

using namespace std;

struct matrix {
    mpz_class m[2][2] = { {1, 1}, {1, 0} };

    matrix friend operator*(const matrix& a, const matrix& b) {
        matrix res = { .m = { {0, 0}, {0, 0} } };

        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                for (int k = 0; k < 2; k++)
                    res.m[i][j] += a.m[i][k] * b.m[k][j];

        return res;
    }

    matrix power(unsigned long long n) {
        if (n == 1)
            return *this;

        matrix ans = { .m = { {1, 0}, {0, 1} } };
        
        while (n) {
            if (n & 1)
                ans = ans * (*this);

            *this = *this * (*this);
            n >>= 1;
        }

        return ans;
    }
};

mpz_class fib(unsigned long long n) {
    if (n <= 2)
        return mpz_class(1);

    auto res = matrix().power(n - 1);
    return res.m[0][0];
}

Conclusion

In my tests, I was able to calculate the 15,000,000th term of the Fibonacci sequence in less than 1 second, anything greater than that will take a long time to calculate, obviously depending on the machine you are using. This obviously is not the best way to calculate the Fibonacci sequence, but it is a good way to learn about linear algebra works and how to use matrices to build fun algorithms.

math
algorithms
linear algebra