Is there a shortcut for taking this matrix to a certain power? One can take a diagonal matrix to a certain power by just taking diagonal elements that power. There is a similar polynomial-time (in respect to the matrix dimensions) shortcut for triangular matrices.
Is there a shortcut for exponentiating a square matrix of the form
$$
\left( \begin{array}{cccc}
a & b & c & d \\
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \end{array} \right)
$$
, where $a$, $b$, $c$, and $d$ are constants?
 A: By using Cayley Hamilton Theorem. If you are doing the calculations by hand, then the problem is simple if the characteristic polynomial has simple roots. 
The general approach is to divide $x^n$ by $x^4-a x^3 - b x^2 - c x - d$ find the remainder. Thus suppose
$$
x^n = Q(x) (x^4-a x^3 - b x^2 - c x - d) + { \alpha x^3 + \beta x^2 + \gamma x + \delta } \tag 1
$$
Here we do not need the quotient $Q(x)$ but just the remainder $\alpha x^3 + \beta x^2 + \gamma x + \delta $
Then, replacing $x$ by $A$ we have
$$
A^n = Q(A) (A^4-a A^3 - b A^2 - c A- d I) + { \alpha A^3 + \beta A^2 + \gamma A + \delta I } \tag 2
$$
In (2) we have $$(A^4-a A^3 - b A^2 - c A- d I) =0 ~\text{Verify this by direct calculation!}$$
so
$$
A^n = { \alpha A^3 + \beta A^2 + \gamma A + \delta I }$$
One can easily calculate
$$A^2=\pmatrix{b+a^2&c+a\,b&d+a\,c&a\,d\cr a&b&c&d\cr 1&0&0&0\cr 0&1&
 0&0\cr }$$
$$A^3=\pmatrix{c+2\,a\,b+a^3&d+a\,c+b^2+a^2\,b&a\,d+b\,c+a^2\,c&
 \left(b+a^2\right)\,d\cr b+a^2&c+a\,b&d+a\,c&a\,d\cr a&b&c&d\cr 1&0&
 0&0\cr }$$
The key to the answer is finding the remainder. One can use recursion (synthetic division) or remainder theorem to find the remainder.
Here is a simple example on finding the remainder when the characteristic equation has simple roots.
Example: Let
$$a=2 , b=1 , c=-2 , d=0  $$
Then the roots of $x^4-a x^3 - b x^2 - c x - d$ are $x=0$, $x=1$, $x=-1$ and $x=-2$.
Now let
$$x^n = Q(x) (x^4-a x^3 - b x^2 - c x - d) + \alpha x^3 + \beta x^2 + \gamma x + \delta$$
We substitute different roots we got as 
Set $x=0$
$$ 0 = \delta \tag 3$$
Set $x=1$:
$$
1 = \alpha+\beta+\gamma+\delta \tag 4$$
Set $x=-1$
$$
(-1)^n = -\alpha+\beta-\gamma+\delta \tag 5$$
Set $x = 2$
$$
2^n = 8 \alpha + 4 \beta + 2 \gamma + \delta \tag6$$
We can solve to get
$$ \alpha={{2^{n}-\left(-1\right)^{n}-3}\over{6}} , \beta={{
 \left(-1\right)^{n}+1}\over{2}} , \gamma=-{{2^{n}+2\,\left(-1\right)
 ^{n}-6}\over{6}} , \delta=0  $$
So we have
$$
A^n = 
{{2^{n}-\left(-1\right)^{n}-3}\over{6}} A^3 +
{{
 \left(-1\right)^{n}+1}\over{2}} A^2 + {{2^{n}+2\,\left(-1\right)
 ^{n}-6}\over{6}} A$$
A: The problem can be solved when the polynomial $x^4 - ax^3 - bx^2 - cx^1 - d$ has 4 distinct (complex) roots.
First we note, that if $\lambda$ is a root, the matrix has eigenvalue $\lambda$ with eigenvector $( \lambda^3, \lambda^2, \lambda, 1 )^T$:
$$
\left( \begin{array}{cccc}
a & b & c & d \\
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \end{array} \right)
\left( \begin{array}{cccc}
\lambda^3  \\
\lambda^2 \\
\lambda \\
1 \end{array} \right)
= 
\left( \begin{array}{cccc}
\lambda^4  \\
\lambda^3 \\
\lambda^2 \\
\lambda \end{array} \right)
= 
\lambda
\left( \begin{array}{cccc}
\lambda^3  \\
\lambda^2 \\
\lambda \\
1 \end{array} \right)
$$
Since it has 4 distinct eigenvalues (and linear independent eigenvectors) we can then diagonalize the matrix.
Let
$$
V = \left( \begin{array}{cccc}
\lambda_1^3 & \lambda_2^3 & \lambda_3^3 & \lambda_4^3  \\
\lambda_1^2 & \lambda_2^2 & \lambda_3^2 & \lambda_4^2\\
\lambda_1 & \lambda_2 & \lambda_3 & \lambda_4\\
1 & 1 & 1 & 1 \end{array} \right)
$$
then we
we have
$$
 D = V^{-1}M V
$$
with the original matrix denoted as $M$ and 
$$
D = \left( \begin{array}{cccc}
\lambda_1 & 0 & 0& 0  \\
0 & \lambda_2 & 0 & 0\\
0& 0 & \lambda_3 & 0\\
0 & 0 & 0 & \lambda_4 \end{array} \right)
$$
Finally, we can compute arbitrary powers of $M$ by
$$
M^n = V D^n V^{-1}
$$
This derivation is analogous to the one of the companion matrix, which is very similar to the present matrix (here different permutation and negative $a,b,c,d$). See http://en.wikipedia.org/wiki/Companion_matrix
