Visual Studio Code – How to calculate nth term of higher order GENERALIZED Fibonacci sequence using matrix exponentiation in C++?

I didn’t know what the code was really trying to accomplish.

So when piecing it together I accidentally rewrote the stuff. Here’s my iterative any-Order fibonacci:

namespace ranges = std::ranges;
namespace views  = ranges::views;
using Int   = uint64_t;

template <size_t Order> static constexpr auto onacci_iterative(uint64_t n) {
    std::array<Int, Order> window{};
    window.back() = 1;
    while (true) {
        if (auto val = window.front())
            if (n-- == 0)
                return val;

        auto sum = ranges::fold_left(window, Int{}, std::plus<>{});
        ranges::copy(views::drop(window, 1), window.begin());
        window.back() = sum;
    }

    throw std::logic_error{"unreachable"};
}

Here’s a simple demo Live Online:

 2: [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987]
 3: [1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504, 927, 1705, 3136, 5768]
 4: [1, 1, 2, 4, 8, 15, 29, 56, 108, 208, 401, 773, 1490, 2872, 5536, 10671]
 5: [1, 1, 2, 4, 8, 16, 31, 61, 120, 236, 464, 912, 1793, 3525, 6930, 13624]
 6: [1, 1, 2, 4, 8, 16, 32, 63, 125, 248, 492, 976, 1936, 3840, 7617, 15109]
 7: [1, 1, 2, 4, 8, 16, 32, 64, 127, 253, 504, 1004, 2000, 3984, 7936, 15808]
 8: [1, 1, 2, 4, 8, 16, 32, 64, 128, 255, 509, 1016, 2028, 4048, 8080, 16128]
 9: [1, 1, 2, 4, 8, 16, 32, 64, 128, 256, 511, 1021, 2040, 4076, 8144, 16272]
10: [1, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1023, 2045, 4088, 8172, 16336]

Matrix Approach

I won’t do all the matrix logic myself, instead using Boost Ublas:

using Matrix = boost::numeric::ublas::matrix<Int>;

First of all, let’s create the exponentation helper like you describe:

template <class T> T my_pow(T b, int64_t e) {
    if (e == 1)
        return b;
    if (e % 2 == 1)
        return prod(b, my_pow(b, e - 1));
    T tmp = my_pow(b, e / 2);
    return prod(tmp, tmp);
}

Next, we need the onacci_matrix generator:

Matrix onacci_matrix(Int order) {
    Matrix m(order, order);
    m.assign(boost::numeric::ublas::zero_matrix<Int>(order, order));

    for (Int i = 0; i < order; ++i)
        m(i, 0) = 1;

    for (Int i = 1; i < order; ++i)
        m(i-1, i) = 1;

    return order % 2 ? trans(m) : m;
}

To check that everything is correct, I created a benchmark comparison:

template <size_t Order> void compare(uint64_t n) {
    auto timed = [](auto f) {
        auto start = now();
        auto r     = f();
        return std::tuple((now() - start), r);
    };

    auto [itime, itval] = timed([n] { return onacci_iterative<Order>(n); });
    auto [mtime, m]     = timed([n] { return my_pow(onacci_matrix(Order), n); });
    auto mval           = m(0, 0);

    fmt::print("order:{} {:<20} {:10L}μs {}n", Order, fmt::format("iterative({}):", n), itime / 1us, itval);
    fmt::print("order:{} {:<20} {:10L}μs {}n", Order, fmt::format("matrix^{}:", n), mtime / 1us, mval);
    fmt::print(" -- itval == mval? {}n", itval == mval);
}

And when you run that like

compare<6>(19); // the sample from the question

It prints

order:6 iterative(19):                0μs 233904
order:6 matrix^19:                   11μs 233904
 -- itval == mval? true

This obviously matches the value from the question. Note that although at n=19 the iterative solution is obviously faster, it’s possible to check validity even without scaling the integer type:

for (Int n = 19; n < 19e9; n *= 10) {
    compare<3>(n);
    compare<4>(n);
    compare<5>(n);
    compare<6>(n);
}

These all correctly check out, 36x "itval == mval? true":

order:3 iterative(19):                0μs 66012
order:3 matrix^19:                    2μs 66012
 -- itval == mval? true
order:4 iterative(19):                0μs 147312
order:4 matrix^19:                    1μs 147312
 -- itval == mval? true
order:5 iterative(19):                0μs 203513
order:5 matrix^19:                    2μs 203513
 -- itval == mval? true
order:6 iterative(19):                0μs 233904
order:6 matrix^19:                    7μs 233904
// ...
order:3 iterative(1900):              4μs 12869305375559491231
order:3 matrix^1900:                  2μs 12869305375559491231
 -- itval == mval? true
order:4 iterative(1900):             13μs 6881832060256714257
order:4 matrix^1900:                  3μs 6881832060256714257
 -- itval == mval? true
order:5 iterative(1900):             13μs 11498825699531939336
order:5 matrix^1900:                  4μs 11498825699531939336
 -- itval == mval? true
order:6 iterative(1900):             14μs 16271437555610019488
order:6 matrix^1900:                  7μs 16271437555610019488
 -- itval == mval? true
// ...
order:3 iterative(1900000000):  3,176,897μs 7052199387128452865
order:3 matrix^1900000000:            5μs 7052199387128452865
 -- itval == mval? true
order:4 iterative(1900000000):  8,792,814μs 84094515927825409
order:4 matrix^1900000000:            7μs 84094515927825409
 -- itval == mval? true
order:5 iterative(1900000000):  8,793,596μs 4174985315637410720
order:5 matrix^1900000000:           10μs 4174985315637410720
 -- itval == mval? true
order:6 iterative(1900000000):  9,794,345μs 7912436391526238456
order:6 matrix^1900000000:           13μs 7912436391526238456
 -- itval == mval? true

Multiprecision

To get the actually correct values, let’s use Boost Multiprecision. Instead of printing the entire values (which will be very large) we print the ₁₀log (basically, the number of digits).

I benchmarked

• Int = uint64_t: (overflows)

  order:6 iterative(190000):        1,364μs 10log:18.4885
  order:6 matrix^190000:                8μs 10log:18.4885
   -- itval == mval? true
  
  real    0m0,014s
  

• Int = mp::int128_t: (overflows)

  order:6 iterative(190000):          803μs 10log:38.4783
  order:6 matrix^190000:               47μs 10log:38.4783
   -- itval == mval? true
  
  real    0m0,011s
  

• Int = mp::int1024_t: (still overflows)

  order:6 iterative(190000):       29,524μs 10log:308.012
  order:6 matrix^190000:              498μs 10log:308.012
   -- itval == mval? true
  
  real    0m0,121s
  

The only types that will likely fit the bill are arbitrary-precision integer types:

• Int = mp::cpp_int: (never overflows, no link dependency)

  order:6 iterative(190000):    1,247,917μs 10log:56515.3
  order:6 matrix^190000:          198,435μs 10log:56515.3
   -- itval == mval? true
  
  real    0m4,039s
  

• Int = mp::mpz_int: (can never overflow, uses GMP/MPFR)

  order:6 iterative(190000):    1,339,042μs 10log:56515.3
  order:6 matrix^190000:           81,846μs 10log:56515.3
   -- itval == mval? true
  
  real    0m4,057s
  

Here’s the code with the mpz_int enabled: Live On Coliru

#define NDEBUG
#include <algorithm>
#include <array>
#include <chrono>
#include <cstdint>
#include <numeric>
#include <ranges>
#include <vector>

using namespace std::chrono_literals;
static constexpr auto now = std::chrono::high_resolution_clock::now;

namespace ranges = std::ranges;
namespace views  = ranges::views;

template <size_t Order, typename Int> static constexpr auto onacci_iterative(uint64_t n) {
    std::array<Int, Order> window{};
    window.back() = 1;
    while (true) {
        if (Int const& val = window.front())
            if (n-- == 0)
                return val;

        auto sum = ranges::fold_left(window, Int{}, std::plus<>{});
        ranges::copy(views::drop(window, 1), window.begin());
        window.back() = sum;
    }

    throw std::logic_error{"unreachable"};
}

#include <fmt/ranges.h>
#include <fmt/ostream.h>

#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/gmp.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>

template <class T> T my_pow(T b, int64_t e) {
    if (e == 1)
        return b;
    if (e % 2 == 1)
        return prod(b, my_pow(b, e - 1));
    T tmp = my_pow(b, e / 2);
    return prod(tmp, tmp);
}

// using Int = uint64_t;
// using Int = boost::multiprecision::int128_t;
// using Int = boost::multiprecision::int1024_t;
// using Int = boost::multiprecision::cpp_int;
using Int    = boost::multiprecision::mpz_int;
using Matrix = boost::numeric::ublas::matrix<Int>;

using Real = boost::multiprecision::mpf_float_50; // just for log10
using boost::multiprecision::log10;

template <> struct fmt::formatter<Int> : fmt::ostream_formatter { };
template <> struct fmt::formatter<Real> : fmt::ostream_formatter { };

Matrix onacci_matrix(size_t order) {
    Matrix m(order, order);
    m.assign(boost::numeric::ublas::zero_matrix<Int>(order, order));

    for (size_t i = 0; i < order; ++i)
        m(i, 0) = 1;

    for (size_t i = 1; i < order; ++i)
        m(i-1, i) = 1;

    return order % 2 ? trans(m) : m;
}

template <size_t Order> void compare(uint64_t n) {
    auto timed = [](auto f) {
        auto start = now();
        auto r     = f();
        return std::tuple((now() - start), r);
    };

    auto [itime, itval] = timed([n] { return onacci_iterative<Order, Int>(n); });
    auto [mtime, m]     = timed([n] { return my_pow(onacci_matrix(Order), n); });
    auto mval           = m(0, 0);

    fmt::print("order:{} {:<20} {:10L}μs 10log:{}n", Order, fmt::format("iterative({}):", n), itime / 1us, log10(Real(itval)));
    fmt::print("order:{} {:<20} {:10L}μs 10log:{}n", Order, fmt::format("matrix^{}:", n), mtime / 1us, log10(Real(mval)));
    fmt::print(" -- itval == mval? {}n", itval == mval);
}

int main() {
    std::locale::global(std::locale("en_US.UTF-8"));

    for (uint64_t n = 19; n < 19e5; n *= 10) {
        compare<3>(n);
        compare<4>(n);
        compare<5>(n);
        compare<6>(n);
    }
}

SUMMARY

So while nothing above addresses your code, it may help you validating your intermediates, and perhaps getting ballpark measures on what performance to expect.