C-Plus-Plus/math/linear_recurrence_matrix.cpp

/**
 * @brief Evaluate recurrence relation using [matrix
 * exponentiation](https://www.hackerearth.com/practice/notes/matrix-exponentiation-1/).
 * @details
 * Given a recurrence relation; evaluate the value of nth term.
 * For e.g., For fibonacci series, recurrence series is `f(n) = f(n-1) + f(n-2)`
 * where `f(0) = 0` and `f(1) = 1`.
 * Note that the method used only demonstrates
 * recurrence relation with one variable (n), unlike `nCr` problem, since it has
 * two (n, r)
 *
 * ### Algorithm
 * This problem can be solved using matrix exponentiation method.
 * @see here for simple [number exponentiation
 * algorithm](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/modular_exponentiation.cpp)
 * or [explaination
 * here](https://en.wikipedia.org/wiki/Exponentiation_by_squaring).
 * @author [Ashish Daulatabad](https://github.com/AshishYUO)
 */
#include <cassert>   /// for assert
#include <cstdint>
#include <iostream>  /// for IO operations
#include <vector>    /// for std::vector STL

/**
 * @namespace math
 * @brief Mathematical algorithms
 */
namespace math {
/**
 * @namespace linear_recurrence_matrix
 * @brief Functions for [Linear Recurrence
 * Matrix](https://www.hackerearth.com/practice/notes/matrix-exponentiation-1/)
 * implementation.
 */
namespace linear_recurrence_matrix {
/**
 * @brief Implementation of matrix multiplication
 * @details Multiplies matrix A and B, given total columns in A are equal to
 * total given rows in column B
 * @tparam T template type for integer as well as floating values, default is
 * long long int
 * @param _mat_a first matrix of size n * m
 * @param _mat_b second matrix of size m * k
 * @returns `_mat_c` resultant matrix of size n * k
 * Complexity: `O(n*m*k)`
 * @note The complexity in this case will be O(n^3) due to the nature of the
 * problem. We'll be multiplying the matrix with itself most of the time.
 */
template <typename T = int64_t>
std::vector<std::vector<T>> matrix_multiplication(
    const std::vector<std::vector<T>>& _mat_a,
    const std::vector<std::vector<T>>& _mat_b, const int64_t mod = 1000000007) {
    // assert that columns in `_mat_a` and rows in `_mat_b` are equal
    assert(_mat_a[0].size() == _mat_b.size());
    std::vector<std::vector<T>> _mat_c(_mat_a.size(),
                                       std::vector<T>(_mat_b[0].size(), 0));
    /**
     * Actual matrix multiplication.
     */
    for (uint32_t i = 0; i < _mat_a.size(); ++i) {
        for (uint32_t j = 0; j < _mat_b[0].size(); ++j) {
            for (uint32_t k = 0; k < _mat_b.size(); ++k) {
                _mat_c[i][j] =
                    (_mat_c[i][j] % mod +
                     (_mat_a[i][k] % mod * _mat_b[k][j] % mod) % mod) %
                    mod;
            }
        }
    }
    return _mat_c;
}
/**
 * @brief Returns whether matrix `mat` is a [zero
 * matrix.](https://en.wikipedia.org/wiki/Zero_matrix)
 * @tparam T template type for integer as well as floating values, default is
 * long long int
 * @param _mat A matrix
 * @returns true if it is a zero matrix else false
 */
template <typename T = int64_t>
bool is_zero_matrix(const std::vector<std::vector<T>>& _mat) {
    for (uint32_t i = 0; i < _mat.size(); ++i) {
        for (uint32_t j = 0; j < _mat[i].size(); ++j) {
            if (_mat[i][j] != 0) {
                return false;
            }
        }
    }
    return true;
}

/**
 * @brief Implementation of Matrix exponentiation
 * @details returns the matrix exponentiation `(B^n)` in `k^3 * O(log2(power))`
 * time, where `k` is the size of matrix (k by k).
 * @tparam T template type for integer as well as floating values, default is
 * long long int
 * @param _mat matrix for exponentiation
 * @param power the exponent value
 * @returns the matrix _mat to the power `power (_mat^power)`
 */
template <typename T = int64_t>
std::vector<std::vector<T>> matrix_exponentiation(
    std::vector<std::vector<T>> _mat, uint64_t power,
    const int64_t mod = 1000000007) {
    /**
     * Initializing answer as identity matrix. For simple binary
     * exponentiation reference, [see
     * here](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/modular_exponentiation.cpp)
     */
    if (is_zero_matrix(_mat)) {
        return _mat;
    }

    std::vector<std::vector<T>> _mat_answer(_mat.size(),
                                            std::vector<T>(_mat.size(), 0));

    for (uint32_t i = 0; i < _mat.size(); ++i) {
        _mat_answer[i][i] = 1;
    }
    // exponentiation algorithm here.
    while (power > 0) {
        if (power & 1) {
            _mat_answer = matrix_multiplication(_mat_answer, _mat, mod);
        }
        power >>= 1;
        _mat = matrix_multiplication(_mat, _mat, mod);
    }

    return _mat_answer;
}

/**
 * @brief Implementation of nth recurrence series.
 * @details Returns the nth term in the recurrence series.
 * Note that the function assumes definition of base cases from `n = 0`
 * (e.g., for fibonacci, `f(0)` has a defined value `0`)
 * @tparam T template type for integer as well as floating values, default is
 * long long int
 * @param _mat [square matrix](https://en.m.wikipedia.org/wiki/Square_matrix)
 * that evaluates the nth term using exponentiation
 * @param _base_cases 2D array of dimension `1*n` containing values which are
 * defined for some n (e.g., for fibonacci, `f(0)` and `f(1)` are defined, and
 * `f(n)` where `n > 1` is evaluated on previous two values)
 * @param nth_term the nth term of recurrence relation
 * @param constant_or_sum_included whether the recurrence relation has a
 * constant value or is evaluating sum of first n terms of the recurrence.
 * @returns the nth term of the recurrence relation in `O(k^3. log(n))`, where k
 * is number of rows and columns in `_mat` and `n` is the value of `nth_term`
 * If constant_or_sum_included is true, returns the sum of first n terms in
 * recurrence series
 */
template <typename T = int64_t>
T get_nth_term_of_recurrence_series(
    const std::vector<std::vector<T>>& _mat,
    const std::vector<std::vector<T>>& _base_cases, uint64_t nth_term,
    bool constant_or_sum_included = false) {
    assert(_mat.size() == _base_cases.back().size());

    /**
     * If nth term is a base case, then return base case directly.
     */

    if (nth_term < _base_cases.back().size() - constant_or_sum_included) {
        return _base_cases.back()[nth_term - constant_or_sum_included];
    } else {
        /**
         * Else evaluate the expression, so multiplying _mat to itself (n -
         * base_cases.length + 1 + constant_or_sum_included) times.
         */
        std::vector<std::vector<T>> _res_matrix =
            matrix_exponentiation(_mat, nth_term - _base_cases.back().size() +
                                            1 + constant_or_sum_included);

        /**
         * After matrix exponentiation, multiply with the base case to evaluate
         * the answer. The answer is always at the end of the array.
         */
        std::vector<std::vector<T>> _res =
            matrix_multiplication(_base_cases, _res_matrix);

        return _res.back().back();
    }
}
}  // namespace linear_recurrence_matrix
}  // namespace math

/**
 * @brief Self test-implementations
 * @returns void
 */
static void test() {
    /*
     * Example 1: [Fibonacci
     * series](https://en.wikipedia.org/wiki/Fibonacci_number);
     *
     * [fn-2    fn-1]  [0      1]  ==   [fn-1   (fn-2 + fn-1)] => [fn-1   fn]
     *                 [1      1]
     *
     * Let A = [fn-2   fn-1], and B = [0   1]
     *                                [1   1],
     *
     * Since, A.B....(n-1 times) = [fn-1   fn]
     * we can multiply B with itself n-1 times to obtain the required value
     */
    std::vector<std::vector<int64_t>> fibonacci_matrix = {{0, 1}, {1, 1}},
                                      fib_base_case = {{0, 1}};

    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               fibonacci_matrix, fib_base_case, 11) == 89LL);
    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               fibonacci_matrix, fib_base_case, 39) == 63245986LL);
    /*
     * Example 2: [Tribonacci series](https://oeis.org/A000073)
     *                    [0   0   1]
     * [fn-3  fn-2  fn-1] [1   0   1]  =  [(fn-2)  (fn-1)  (fn-3 + fn-2 + fn-1)]
     *                    [0   1   1]
     *                                 => [fn-2     fn-1    fn]
     *
     *                                       [0   0   1]
     * Let A = [fn-3   fn-2   fn-1], and B = [1   0   1]
     *                                       [0   1   1]
     *
     * Since, A.B....(n-2 times) = [fn-2  fn-1   fn]
     * we will have multiply B with itself n-2 times to obtain the required
     * value ()
     */

    std::vector<std::vector<int64_t>> tribonacci = {{0, 0, 1},
                                                    {1, 0, 1},
                                                    {0, 1, 1}},
                                      trib_base_case = {
                                          {0, 0, 1}};  // f0 = 0, f1 = 0, f2 = 1

    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               tribonacci, trib_base_case, 11) == 149LL);
    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               tribonacci, trib_base_case, 36) == 615693474LL);

    /*
     * Example 3: [Pell numbers](https://oeis.org/A000129)
     * `f(n)  = 2* f(n-1) + f(n-2); f(0) = f(1) = 2`
     *
     * [fn-2  fn-1] [0   1]  =  [(fn-1)  fn-2 + 2*fn-1)]
     *              [1   2]
     *                       => [fn-1     fn]
     *
     * Let A = [fn-2  fn-1], and B = [0   1]
     *                               [1   2]
     */

    std::vector<std::vector<int64_t>> pell_recurrence = {{0, 1}, {1, 2}},
                                      pell_base_case = {
                                          {2, 2}};  // `f0 = 2, f1 = 2`

    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               pell_recurrence, pell_base_case, 15) == 551614LL);
    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               pell_recurrence, pell_base_case, 23) == 636562078LL);

    /*
     * Example 4: Custom recurrence relation:
     * Now the recurrence is of the form `a*f(n-1) + b*(fn-2) + ... + c`
     * where `c` is the constant
     * `f(n)  = 2* f(n-1) + f(n-2) + 7; f(0) = f(1) = 2, c = 7`
     *
     *                   [1   0   1]
     * [7,  fn-2,  fn-1] [0   0   1]
     *                   [0   1   2]
     * =  [7,  (fn-1),  fn-2 + 2*fn-1) + 7]
     *
     * => [7,    fn-1,     fn]
     * :: Series will be 2, 2, 13, 35, 90, 222, 541, 1311, 3170, 7658, 18493,
     * 44651, 107802, 260262, 628333, 1516935, 362210, 8841362, 21344941,
     * 51531251
     *
     * Let A = [7,  fn-2,  fn-1], and B = [1   0   1]
     *                                    [0   0   1]
     *                                    [0   1   2]
     */

    std::vector<std::vector<int64_t>>
        custom_recurrence = {{1, 0, 1}, {0, 0, 1}, {0, 1, 2}},
        custom_base_case = {{7, 2, 2}};  // `c = 7, f0 = 2, f1 = 2`

    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               custom_recurrence, custom_base_case, 10, 1) == 18493LL);
    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               custom_recurrence, custom_base_case, 19, 1) == 51531251LL);

    /*
     * Example 5: Sum fibonacci sequence
     * The following matrix evaluates the sum of first n fibonacci terms in
     * O(27. log2(n)) time.
     * `f(n) = f(n-1) + f(n-2); f(0) = 0, f(1) = 1`
     *
     *                           [1   0   0]
     * [s(f, n-1),  fn-2,  fn-1] [1   0   1]
     *                           [1   1   1]
     *   => [(s(f, n-1)+f(n-2)+f(n-1)), (fn-1),  f(n-2)+f(n-1)]
     *
     *   => [s(f, n-1)+f(n),    fn-1,     fn]
     *
     *   => [s(f, n),    fn-1,     fn]
     *
     * Sum of first 20 fibonacci series:
     * 0, 1, 2, 4, 7, 12, 20, 33, 54, 88, 143, 232, 376, 609, 986, 1596, 2583,
     * 4180, 6764
     *          f0  f1  s(f,1)
     * Let A = [0    1    1], and B = [0   1   1]
     *                                [1   1   1]
     *                                [0   0   1]
     */

    std::vector<std::vector<int64_t>> sum_fibo_recurrence = {{0, 1, 1},
                                                             {1, 1, 1},
                                                             {0, 0, 1}},
                                      sum_fibo_base_case = {
                                          {0, 1, 1}};  // `f0 = 0, f1 = 1`

    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               sum_fibo_recurrence, sum_fibo_base_case, 13, 1) == 609LL);
    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               sum_fibo_recurrence, sum_fibo_base_case, 16, 1) == 2583LL);
    /*
     * Example 6: [Tribonacci sum series](https://oeis.org/A000073)
     *                               [0   0   1   1]
     * [fn-3  fn-2  fn-1  s(f, n-1)] [1   0   1   1]
     *                               [0   1   1   1]
     *                               [0   0   0   1]
     *
     * = [fn-2, fn-1, fn-3 + fn-2 + fn-1, (fn-3 + fn-2 + fn-1 + s(f, n-1))]
     *
     * => [fn-2, fn-1, fn,  fn + s(f, n-1)]
     *
     * => [fn-2, fn-1, fn, s(f, n)]
     *
     * Sum of the series is: 0, 0, 1, 2, 4, 8, 15, 28, 52, 96, 177, 326, 600,
     * 1104, 2031, 3736, 6872, 12640, 23249, 42762
     *
     * Let A = [fn-3   fn-2   fn-1   s(f, n-1)], and
     *     [0   0   1   1]
     * B = [1   0   1   1]
     *     [0   1   1   1]
     *     [0   0   0   1]
     *
     * Since, A.B....(n-2 times) = [fn-2  fn-1   fn]
     * we will have multiply B with itself n-2 times to obtain the required
     * value
     */

    std::vector<std::vector<int64_t>> tribonacci_sum = {{0, 0, 1, 1},
                                                        {1, 0, 1, 1},
                                                        {0, 1, 1, 1},
                                                        {0, 0, 0, 1}},
                                      trib_sum_base_case = {{0, 0, 1, 1}};
    // `f0 = 0, f1 = 0, f2 = 1, s = 1`

    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               tribonacci_sum, trib_sum_base_case, 18, 1) == 23249LL);
    assert(math::linear_recurrence_matrix::get_nth_term_of_recurrence_series(
               tribonacci_sum, trib_sum_base_case, 19, 1) == 42762LL);
}

/**
 * @brief Main function
 * @returns 0 on exit
 */
int main() {
    test();  // run self-test implementations
    return 0;
}