In this guide, I'd like to share with you a technique I use to solve a wide variety of problems that requires matrix exponentiation. The steps to use this technique will be provided along with an example of usage in a variation of the problem TETRAHRD — Sum of Tetranacci numbers from SPOJ. In this example, let's consider we want the sum of the interval of the Tetranacci numbers from 0 to a number N, given in the input.

Steps

Create a DP solution

Firstly, use dynamic programming to solve it. Use two vectors, one representing the current state and the other representing the next state. The next_state vector must always preserve the following structure:

next_state[idx] = c_1 * cur_state[idx_1] + c_2 * cur_state[idx_2] + ... + c_k * cur_state[idx_k]

For this example, we'll use a trick to compute the answer. Instead of creating a variable and summing it every time to the current value of the Tetranacci number, we'll use one index from the DP array to store that.

// t(i) = t(i - 1) + t(i - 2) + t(i - 3) + t(i - 4)
// t(0) = t(1) = t(2) = 0, t(3) = 1

vector<int> cur_state(5), next_state(5);

// cur_state from 0 to 3, represents the Tetranacci number of n - 3, n - 2, n -
// 1, n. Initially they represent t(0), t(1), t(2), t(3). Index 4 represents
// the prefix sum up to n - 1. Initially, it's equal to t(0) + t(1) + t(2) = 0.
cur_state[3] = 1;

for (int i = 4; i <= n; ++i) {
  // t(i) = t(i - 1) + t(i - 2) + t(i - 3) + t(i - 4)
  next_state[3] =
      1 * cur_state[3] + 1 * cur_state[2] + 1 * cur_state[1] + 1 * cur_state[0];
  // next_state[4] = t(i - 1) + t(i - 2) + ... + t(1) + t(0)
  next_state[4] = 1 * cur_state[4] + 1 * cur_state[3];

  for (int j = 0; j < 3; ++j)
    // shifting the values of the Tetranacci number to the left to be used in
    // the next iteration.
    next_state[j] = 1 * cur_state[j + 1];

  cur_state = next_state;
}

Create the matrix

Just repeat the coefficients according to the DP solution.

vector<vector<int>> mat(5, vector<int>(5, 0));

// t(i) = 1 * t(i - 1) + 1 * t(i - 2) + 1 * t(i - 3) + 1 * t(i - 4) 
mat[3][3] = 1, mat[3][2] = 1, mat[3][1] = 1, mat[3][0] = 1;
// ans += t(i - 1)
mat[4][4] = 1, mat[4][3] = 1;

// next_state[j] = 1 * cur_state[j + 1];
for (int i = 0; i < 3; ++i)
  mat[i][i + 1] = 1;

Exponentiate the matrix

We need to exponentiate the matrix according to the number of iterations in the calculation of the DP. In this case (n - k + 1). k represents the number of coefficients in the expression (in this case, it's equal to 4). We need to add one since the (n + 1)-th state will contain the answer for the n-th state.

mat = mat.exp(n - k + 1); // let's assume we have this method for the matrix
                          // which exponentiates it to the power of (n - k + 1).

Compute the answer

Since we already have the final values for the coefficients, we only need to compute them based on the first known state.

vector<vector<int>> ini(5, vector<int>(1, 0));
// only the value of a[3] will be equal to 1
ini[3][0] = 1;

mat = mat * ini;

In the end, the answer will be at the same index as the DP approach.

cout << "The answer is " << mat[4][0] << endl;

My code to solve TETRAHRD

Problems

Jzzhu and Sequences
SEQ — Recursive Sequence
Tetrahedron
TETRAHRD — Sum of Tetranacci numbers
Wet Shark and Blocks
FLIB — Flibonakki