If you are not interested in linear recurrences, or are already aware of Cayley-Hamilton theorem, you can probably stop reading now.

Suppose we want to solve the following recurrence:

(T[0], T[1], T[2]) = (1, 2, 3), and

T[n] = T[n - 1] + 3 T[n - 2] + 8 T[n - 3], for n >= 3.

The usual way to solve such recurrence is via matrix exponentiation. We create a matrix M

1 3 8 M = 1 0 0 0 1 0 v = 3 2 1

In order to find T[n] we compute M^{n - 2}, and take the dot product of the first row of M with the vector v. This gives us an O (k^3 lg n) algorithm assuming that the size of the matrix was (k x k), and we are using the most elementary algorithm for matrix multiplication.

**Cayley-Hamilton Theorem:** Any square matrix is a solution of its characteristic polynomial. In our example, the characteristic polynomial is:

| (x - 1) -3 -8 | | -1 x 0 | | 0 -1 x | = x^3 - x^2 - 3x - 8

Note that, the characteristic polynomial looks very similar to the recurrence we want to solve. This is no coincidence, it can be shown that the characteristic polynomial of the matrix M corresponding to a linear recurrence looks the same as the original recurrence. So we do not need to calculate the symbolic determinant of the matrix, but we can use the recurrence to compute this polynomial.

Now, according to Cayley-Hamilton theorem:

M^{3} - M^{2} - 3M - 8I = 0, i.e.,

M^{3} = M^{2} + 3M + 8I

Here I is the identity matrix.

Let’s calculate some higher powers of M.

M^{4} = M (M^{2} + 3M + 8I) = M^{3} + 3M^{2} + 8M

= (M^{2} + 3M + 8I) + 3M^{2} + 8M = 4M^{2} + 11M + 8I

Similarly,

M^{5} = 15M^{2} + 20M + 32I

One can observe that for any n, the matrix M^{n} can be represented as a linear combination of M^{2}, M and I. Also, if we know this representation for M^{n}, we can easily calculate the representation for M^{n+1}.

M^{n} = aM^{2} + bM + cI

==> M^{n + 1} = M (aM^{2} + bM + cI) = (a + b) M^{2} + (3a + c) M + 8aI

For a (k x k) matrix, the characteristic polynomial will be of degree k, and the above step can be performed in O (k) time, i.e., we can compute the representation of M^{n+1} from M^{n} in O (k) time.

Now, let us see if we can calculate the representation of M^{2n} from the representation of M^{n}.

M^{2n} = (aM^{2} + bM + cI ) (aM^{2} + bM + cI )

= a^{2} M^{4} + 2ab M^{3} + (b^{2} + 2ac) M^{2} + 2bc M + c^{2} I

We can replace M^{4} and M^{3} by their corresponding representation that we have already calculated. This means, in order to compute M^{2n} from M^{n}, we need to multiply two degree k polynomials, and then replace M^{2k -2}, M^{2k - 3}, …, M^{k} by their representation. This can be done in O (k^2) time using elementary methods.

Now, we know how to compute the representation of:

- M
^{n+1}from M^{n}, and - M
^{2n}from M^{n}

Both steps have O (k^2) complexity. Since these are the only steps needed for binary exponentiation algorithm, we can compute the representation of M^{n} in O (k^2 lg n) time.

Now, to compute the actual matrix M^{n}, we still need to add these k matrices which are present in its representation. However, if you notice, we do not need the whole matrix M^{n}, but only its first row, which again can be computed in O (k^2) time.

Computing the first row of M^{2}, M^{3}, …, M^{k - 1} is not difficult either, and can be done is O (k^2) time.

This gives us an O (k^2 lg n) algorithm to solve the linear recurrence, which is faster by a factor k, compared to matrix exponentiation method. The difference probably will not make much difference if k is small. However, for large matrices the difference is significant, e.g., project euler 258