Calculate huge Fibonacci Numbers in C++ using GMP and boost::multiprecision

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.

• 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 – Juho Feb 8 at 22:58