Matrix Exponentiation in C++

In this tutorial, we will learn about Matrix Exponentiation and apply it to solve linear recurrences in one variable with C++.

Given a linear recurrence relation in one variable along with the base cases, the task is to find the Nth term in logarithmic time complexity. Since the Nth term can be very large, we will compute it modulo 109 + 7 (a big prime number).

Throughout this tutorial, we will work on the following recurrence relation: –

            RN = 2 * RN – 1 + RN – 2 + 3 * RN – 3

            [R0 = 1, R1 = 2, R2 = 3]

Examples

  1. Input: N = 3
    Output: 11
    Explanation:    R3 = 2 * R2 + R1 + 3 * R0 = 2 * 3 + 2 + 3 * 1 = 11.  

  2. Input: N = 4
    Output: 31
    Explanation:    R4 = 2 * R3 + R2 + 3 * R1 = 2 * 11 + 3 + 3 * 2 = 31.

  3. Input: N = 5
    Output: 82
    Explanation:    R5 = 2 * R4 + R3 + 3 * R2 = 2 * 31 + 11 + 3 * 3 = 82.

Terminologies

Recurrence Relation: A recurrence relation is an algebraic equation that defines a mathematical sequence where the next term of the sequence is a polynomial function of its previous terms.

Linear Recurrence Relation: It is a recurrence relation where the sequence is defined by a degree 1 polynomial. In other words, the next term of the sequence is a linear function of its previous terms. The problems on linear recurrence relations in one variable can be easily solved using Matrix Exponentiation.

Degree of a Linear Recurrence Relation: It is the no. of previous terms that define the next term of the sequence. 

Matrix Exponentiation: It is the process of exponentiating a matrix of size k x k to the power N in O (k3 log N) time complexity.

Coefficient Matrix: It is the matrix that describes a linear recurrence relation in one variable.

Maths Behind The Algorithm

The given linear recurrence relation can be written as: –

Matrix Exponentiation

Here, C is known as the coefficient matrix. It completely depends on the linear recurrence relation.

The coefficient matrix C can be exponentiated to the power N – 2 in O (log N). The approach is pretty similar to modular exponentiation. Refer this for more details.

Refer to this for more details on Matrix Multiplication.

Now, let A = CN – 2. Then: –

RN = 3 * A00 + 2 * A01 + 1 * A02

Algorithm

  • Convert the given linear recurrence relation to matrix form by defining the coefficient matrix.

  • Exponentiate the coefficient matrix to the power N – 2 using Matrix Exponentiation.

  • Find the Nth term by using the exponentiated coefficient matrix and the base cases.

C++ Implementation of the above Algorithm of Matrix Exponentiation

#include <bits/stdc++.h>
using namespace std;

// No of terms in the Recurrence Relation.
const int N = 3;

const long long M = 1000000007;

// Multiplies two matrices A and B and stores the result in A.
void multiply (long long A[N][N], long long B[N][N])
{
    long long R[N][N];
    
    // Multiply A and B and store result in R.
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            R[i][j] = 0;
            
            for (int k = 0; k < N; k++)
            {
                R[i][j] = (R[i][j] + A[i][k] * B[k][j]) % M;
            }
        }
    }
    
    // Copy contents of R in A.
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            A[i][j] = R[i][j];
        }
    }
}

// Raise matrix A to the power of n in O(log n).
void power_matrix (long long A[N][N], int n)
{
    long long B[N][N];
    
    // B = Identity Matrix.
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            B[i][j] = A[i][j];
        }
    }
    
    // A = A * A ^ (n - 1).
    n = n - 1;
    while (n > 0)
    {
        // If n is odd, A = A * B.
        if (n & 1)
            multiply (A, B);
            
        // B = B * B.
        multiply (B,B);
        
        // n = n / 2.
        n = n >> 1;  
    }
}

// A = Coefficient Matrix, B = Base Matrix. 
// It returns the nth term of the recurrence relation formed from A and B in O(log n).
long long solve_recurrence (long long A[N][N], long long B[N][1], int n)
{
    //Base Cases.
    if (n < N)
        return B[N - 1 - n][0];
    
    // A = A ^ (n - N + 1).
    power_matrix (A, n - N + 1);
    
    long long result = 0;
    
    for (int i = 0; i < N; i++)
        result = (result + A[0][i] * B[i][0]) % M;
    
    return result;
}

// Driver Code.
int main ()
{
    /* 
        The recurrence relation used here is: -
        R(n) = 2 * R(n-1) + R(n-2) + 3 * R(n-3).
        Base Cases: R(0) = 1, R(1) = 2, R(2) = 3.
    */
    
    // Forming the Coefficient Matrix
    long long A[N][N] = {{2, 1, 3}, {1, 0, 0}, {0, 1, 0}};
    
    //Forming the Base Matrix
    long long B[N][1] = {{3}, {2}, {1}};
    
    int n = 5;
    
    long long R_n = solve_recurrence (A, B, n);
    
    cout << "R_" << n << " = " << R_n; 

    return 0;
}

Output

R_5 = 82

Complexity Analysis

Matrix multiplication takes O (k3) when the matrix is of size k x k. In our case, k is 3. And since 33 is expected to be very small compared to N, we can say that the matrices are getting multiplied in O (1). The matrix as a whole can be exponentiated in O (log N).
Therefore, the time complexity of the above code is O (log N), although the time complexity of the algorithm is O (k3 log N), where k is the degree of the linear recurrence relation in one variable.

The coefficient matrix in our case is a 3 x 3 matrix, which is independent of N.
Therefore, the space complexity of the above code is O (1).

Leave a Reply

Your email address will not be published.