3
\$\begingroup\$

Problem Statement

Calculate the exact value of the n-th Fibonacci number, say the one-billionth.

Algorithm

The algorithm is based on the idea that Fibonacci numbers can be represented as 2x2 matrices, in particular as powers of the matrix [ 1 1; 1 0]. Matrices of this form can be represented implicitly without having to store all 4 elements as [ a+b a; a b ]. In particular, we only need to store a and b to represent the matrix. This not only saves memory but we can perform fewer operations when multiplying or squaring matrices if we work on the implicit form.

The other aspect of the algorithm is exponentiation by squaring, the idea that we can calculate a^(2^n) very quickly by starting with a and repeatedly squaring. To calculate the n-th Fibonacci number, we simply compute the n-th power of the basis matrix (using the implicit form) and read off the diagonal element a.

Build and Implementation Details

For performance we use the boost::multiprecision::mpz_int wrapper around the Gnu MP library. The mpz_int type is a lightweight wrapper around GMP which gives it a C++-style interface with overloaded operators. The GMP is an "arbitrary precision arithmetic" library for doing math with numbers that don't fit into 64 bits. I installed these libraries on Debian/Ubuntu like so:

sudo apt install libboost-all-dev libgmp-dev

The following program was built like so:

g++ -std=c++17 -O3 -o fib main.cpp -lgmp

Program

#include <iostream>
#include <boost/multiprecision/gmp.hpp>

typedef boost::multiprecision::mpz_int bigint;

namespace fib {

class ImplicitMatrix {
public:
    bigint a;
    bigint b;

    ImplicitMatrix(const ImplicitMatrix& other)
        : a(other.a)
        , b(other.b)
    {
    }

    ImplicitMatrix(
        const bigint& _a,
        const bigint& _b)
        : a(_a)
        , b(_b)
    {
    }
};

ImplicitMatrix multiply(
    const ImplicitMatrix& x,
    const ImplicitMatrix& y)
{
    return {
        y.a * (x.a + x.b) + x.a * y.b,
        x.a * y.a + x.b * y.b
    };
}

ImplicitMatrix square(const ImplicitMatrix& x)
{
    // we save one bigint multipication by
    // specializing for the case of squaring.
    bigint a2 = x.a * x.a;
    bigint _2ab = (x.a * x.b) << 1;
    bigint b2 = x.b * x.b;
    return {
        a2 + _2ab,
        a2 + b2
    };
}

ImplicitMatrix power(
    const ImplicitMatrix& x,
    const bigint& m)
{
    if (m == 0) {
        return { 0, 1 };
    }
    else if (m == 1) {
        return x;
    }

    // powers of two by iterated squaring
    ImplicitMatrix powerOfTwo = x;
    bigint n = 2;
    while (n <= m) {
        powerOfTwo = square(powerOfTwo);
        n = n * 2;
    }
    // recurse for remainder
    ImplicitMatrix remainder = power(x, m - n / 2);

    return multiply(powerOfTwo, remainder);
}

} // end namespace fib

bigint fibonacci(const bigint& n)
{
    fib::ImplicitMatrix f1{ 1, 0 };
    return fib::power(f1, n).a;
}

int main(int argc, const char* argv[])
{
    // require the first command line argument
    if (argc < 2) {
        std::cerr << "USAGE: fib N" << std::endl;
        return 1;
    }
    // parse the first argument into a bigint
    bigint n(argv[1]);

    // force the calculation of the n-th fibonacci
    // number but do not print it.
    volatile bigint f = fibonacci(n);
    return 0;
}

I am looking for style advice, modern C++ tips, performance suggestions, and best practices.

\$\endgroup\$
  • \$\begingroup\$ I don't see the point of defining your own class here, it's just a placeholder for two (public) member variables. I'd rather just typedef your ImplicitMatrix as say std::pair \$\endgroup\$ – Juho Feb 8 at 22:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.