# COME ON CODE ON

A blog about programming and more programming.

## Recurrence Relation and Matrix Exponentiation

Recurrence relations appear many times in computer science. Using recurrence relation and dynamic programming we can calculate the nth term in O(n) time. But many times we need to calculate the nth in O(log n) time. This is where Matrix Exponentiation comes to rescue.

We will specifically look at linear recurrences. A linear recurrence is a sequence of vectors defined by the equation Xi+1 = M Xi for some constant matrix M. So our aim is to find this constant matrix M, for a given recurrence relation.

Let’s first start by looking at the common structure of our three matrices Xi+1, Xi and M. For a recurrence relation where the next term is dependent on last k terms, Xi+1 and Xi+1 are matrices of size 1 x k and M is a matrix of size k x k.

```| f(n+1)   |       | f(n)   |
|  f(n)    |       | f(n-1) |
| f(n-1)   | = M x | f(n-2) |
| ......   |       | ...... |
| f(n-k+1) |       | f(n-k) |```

Let’s look at different type of recurrence relations and how to find M.

1. Let’s start with the most common recurrence relation in computer science, The Fibonacci Sequence. Fi+1 = Fi + Fi-1.

```| f(n+1) | = M x | f(n)   |
| f(n)   |       | f(n-1) |```

Now we know that M is a 2 x 2 matrix. Let it be

```| f(n+1) | = | a b | x | f(n)   |
| f(n)   |   | c d |   | f(n-1) |```

Now a*f(n) + b*f(n-1) = f(n+1) and c*f(n) + d*f(n-1) = f(n). Solving these two equations we get a=1,b=1,c=1 and d=0. So,

```| f(n+1) | = | 1 1 | x | f(n)   |
| f(n)   |   | 1 0 |   | f(n-1) |```

2. For recurrence relation f(n) = a*f(n-1) + b*f(n-2) + c*f(n-3), we get

```| f(n+1) | = | a b c | x | f(n)   |
| f(n)   |   | 1 0 0 |   | f(n-1) |
| f(n-1) |   | 0 1 0 |   | f(n-2) |```

3. What if the recurrence relation is f(n) = a*f(n-1) + b*f(n-2) + c, where c is a constant. We can also add it in the matrices as a state.

```| f(n+1) | = | a b 1 | x | f(n)   |
| f(n)   |   | 1 0 0 |   | f(n-1) |
| c      |   | 0 0 1 |   | c      |```

4. If a recurrence relation is given like this f(n) = f(n-1) if n is odd, f(n-2) otherwise, we can convert it to f(n) = (n&1) * f(n-1) + (!(n&1)) * f(n-2) and substitute the value accordingly in the matrix.

5.If there are more than one recurrence relation, g(n) = a*g(n-1) + b*g(n-2) + c*f(n), and, f(n) = d*f(n-1) + e*f(n-2). We can still define the matrix X in following way

```| g(n+1) |   | a b c 0 |   | g(n)   |
| g(n)   | = | 1 0 0 0 | x | g(n-1) |
| f(n+2) |   | 0 0 d e |   | f(n+1) |
| f(n+1) |   | 0 0 1 0 |   | f(n)   |```

Now that we have got our matrix M, how are we going to find the nth term.
Xi+1 = M Xi
(Multiplying M both sides)
M * Xi+1 = M * M Xi
Xi+2 = M^2 Xi
..
Xi+k = M^k Xi

So all we need now is to find the matrix M^k to find the k-th term. For example in the case of Fibonacci Sequence,

```M^2 = | 2 1 |
| 1 1 |```

Hence F(2) = 2.

Now we need to learn to find M^k in O(n3 log b) time. The brute force approach to calculate a^b takes O(b) time, but using a recursive divide-and-conquer algorithm takes only O(log b) time:

• If b = 0, then the answer is 1.
• If b = 2k is even, then ab = (ak)2.
• If b is odd, then ab = a * ab-1.

We take a similar approach for Matrix Exponentiation. The multiplication part takes the O(n3) time and hence the overall complexity is O(n3 log b). Here’s a sample code in C++ using template class:

```#include<iostream>
using namespace std;

template< class T >
class Matrix
{
public:
int m,n;
T *data;

Matrix( int m, int n );
Matrix( const Matrix< T > &matrix );

const Matrix< T > &operator=( const Matrix< T > &A );
const Matrix< T > operator*( const Matrix< T > &A );
const Matrix< T > operator^( int P );

~Matrix();
};

template< class T >
Matrix< T >::Matrix( int m, int n )
{
this->m = m;
this->n = n;
data = new T[m*n];
}

template< class T >
Matrix< T >::Matrix( const Matrix< T > &A )
{
this->m = A.m;
this->n = A.n;
data = new T[m*n];
for( int i = 0; i < m * n; i++ )
data[i] = A.data[i];
}

template< class T >
Matrix< T >::~Matrix()
{
delete [] data;
}

template< class T >
const Matrix< T > &Matrix< T >::operator=( const Matrix< T > &A )
{
if( &A != this )
{
delete [] data;
m = A.m;
n = A.n;
data = new T[m*n];
for( int i = 0; i < m * n; i++ )
data[i] = A.data[i];
}
return *this;
}

template< class T >
const Matrix< T > Matrix< T >::operator*( const Matrix< T > &A )
{
Matrix C( m, A.n );
for( int i = 0; i < m; ++i )
for( int j = 0; j < A.n; ++j )
{
C.data[i*C.n+j]=0;
for( int k = 0; k < n; ++k )
C.data[i*C.n+j] = C.data[i*C.n+j] + data[i*n+k]*A.data[k*A.n+j];
}
return C;
}

template< class T >
const Matrix< T > Matrix< T >::operator^( int P )
{
if( P == 1 ) return (*this);
if( P & 1 ) return (*this) * ((*this) ^ (P-1));
Matrix B = (*this) ^ (P/2);
return B*B;
}

int main()
{
Matrix<int> M(2,2);
M.data[0] = 1;M.data[1] = 1;
M.data[2] = 1;M.data[3] = 0;

int F[2]={0,1};
int N;
while (~scanf("%d",&N))
if (N>1)
printf("%lld\n",(M^N).data[0]);
else
printf("%d\n",F[N]);
}```

Written by fR0DDY

May 8, 2011 at 5:26 PM

Posted in Programming

### 14 Responses

1. Nice job…well explained

Ranjan

May 18, 2011 at 12:53 PM

2. Good explanation. I have just learnt a new thing. An interesting one, indeed.

Just a question. Can not the matrix exponentiation be done with the diagonalization method? It would be faster, although I am not very sure. I am a little sleepy right now.

Areku

May 22, 2011 at 3:53 AM

3. Just a question. Can not the matrix exponentiation be done with the diagonalization method? It would be faster, although I am not very sure. I am a little sleepy right now.
+1

Jabez

June 5, 2011 at 1:43 AM

4. hi how your code is high lighted in the wordpress
can you please tell me up??
i am not getting plugin option also in my word press did you paid for this??

gauravalgo

June 17, 2011 at 5:33 PM

5. […] Using Matrix Multiplication In the last post we learned how to use Fast Matrix Multiplication to calculate functions having linear equations in […]

6. i want to know that can Determinant of matrix be find out without using any floating point operations.

gauravalgo

August 16, 2011 at 3:59 PM

7. your code in pyton and simple c is better in understanding(also many coders understand like me) than c++ template so please make your write ups in these languages

gauravalgo

August 27, 2011 at 10:09 PM

“2. For recurrence relation f(n) = a*f(n-1) + b*f(n-2) + c*f(n-3)”
how do we calculate a,b,c by applying (subset sum algorithm) or by some else how?

and second thing that your code in c++ template so may be
cannot undestand by many users so i have this techical writing for those

http://algosolver.blogspot.com/2011/09/fibonacci-log-n.html

enjoy ; )

gauravalgo

September 1, 2011 at 6:00 PM

9. The first K terms of the sequence F(0), F(1), F(2), … , F(K-1) are given and the recurrence is, F(n) = F(n-1) + F(n-2) + … + F(n-K) for n >= K. Given N, find the value of F(N).

vcfvamcvmc

October 27, 2011 at 2:29 AM

10. @fRODDY….i din understand how u’ve written x_i+1 matrix as g_n+1,g_n,f_n+2,f_n+1 ???why isn’t it f_n+1,f_n????please clarify???how can we deduce the no of elements in the k*1 matrix if reccurences like this are given????

Lalit

May 16, 2012 at 2:02 PM

11. Can we calculate this recurrence equation using matrix ????
f(n)=f(n-1)+f(n-2)-2^n

Nafis Islam

September 10, 2012 at 1:40 AM

12. Can anyone tell me how the matrix turns out to be {2,2,1,0} for f(n) = f(n-1)+4f(n-2)+2f(n-3) in the below link?
http://www.codechef.com/viewsolution/547362

unknown

July 15, 2013 at 10:52 PM

13. can u give me the exact order of matrices! ie the matrice (Xi+1) and matric(Xi) ‘s order to be use while programming this approach

paralized

September 27, 2014 at 2:10 PM