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 N^th term in logarithmic time complexity. Since the N^thterm can be very large, we will compute it modulo 10⁹ + 7 (a big prime number).

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

R_N = 2 * R_{N – 1} + R_{N – 2} + 3 * R_{N – 3}

[R₀ = 1, R₁ = 2, R₂ = 3]

Examples

Input: N = 3
Output: 11
Explanation: R₃ = 2 * R₂ + R₁ + 3 * R₀ = 2 * 3 + 2 + 3 * 1 = 11.
Input: N = 4
Output: 31
Explanation: R4 = 2 * R₃ + R₂ + 3 * R₁ = 2 * 11 + 3 + 3 * 2 = 31.
Input: N = 5
Output: 82
Explanation: R5 = 2 * R₄ + R₃ + 3 * R₂ = 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 (k³ 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: –

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 = C^{N – 2}. Then: –

R_N = 3 * A₀₀ + 2 * A₀₁ + 1 * A₀₂

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 N^th 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 (k³) when the matrix is of size k x k. In our case, k is 3. And since 3³ 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 (k³ 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).