Computing High Powers of a Matrix With Polynomial Entries I have a matrix where the terms are from a polynomial with two variables. Specifically it's
$C = \left(\begin{array}{cccc}
1 & 0 & 1 & 0 \\
s &  0 & s & 0 \\
0 & 1 & 0 & 1 \\
0 & st & 0 & st 
\end{array}\right)$ 
Is there a good method for computing (exactly) the n-th power of this matrix?
 A: In general, the best way to calculate arbitrary $n$-th powers of matrices is to diagonalize them, if possible. Writing $C = SDS^{-1}$, then $C^n = SD^nS^{-1}$.
Your matrix can be diagonalized. The calculation can be made a little less gruesome by scaling the 2nd and 4th ordinates by a factor of $1/s$, here.
A: As noted by @Simon S, the matrix can be diagonalized. By the spectral theorem
$$
%
 \mathbf{A} = \mathbf{P}^{-1} \mathbf{\Lambda} \, \mathbf{P} \\
%
$$
Instead of multiplying with the full matrix $\mathbf{A}$, work with the diagonal matrix of eigenvalues. For $k\in\mathbb{N}$:
$$
 \mathbf{A}^{k} = \mathbf{P}^{-1} \mathbf{\Lambda}^{k} \, \mathbf{P} \\
$$

The eigenvalue spectrum is
$$
  \lambda \left( \mathbf{A} \right) = \left\{
\frac{1}{2} \left(1 + st \pm \sqrt{s^2 t^2-2 s t+4 s+1} \right), 0, 0
\right\}
$$
The matrix of eigenvalues
$$
\Lambda =
\left(
\begin{array}{cccc}
 \lambda_{+} & 0 & 0 & 0 \\
 0 & \lambda_{-} & 0 & 0 \\
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
\end{array}
\right)
$$
The matrix of eigenvectors
$$
\begin{align}
%
\mathbf{P} &=
\left(
\begin{array}{cccc}
 -\frac{\lambda_{+}}{2 s^2 t} & -\frac{\lambda_{-}}{2 s^2 t} & -1 & 0 \\
 -\frac{\lambda_{+}}{2 s t} & -\frac{\lambda_{-}}{2 s t} & 0 & -1 \\
 \frac{1}{s t} & \frac{1}{s t} & 1 & 0 \\
 1 & 1 & 0 & 1 \\
\end{array}
\right), \\[4pt]
%
\mathbf{P}^{-1} &=
\left(
\begin{array}{cccc}
 \frac{s t (2 s t-\lambda_{-})}{(t-1) (\lambda_{-}-\lambda_{+})} & \frac{t (\lambda_{-}-2 s)}{(t-1) (\lambda_{-}-\lambda_{+})} & \frac{s t (2 s t-\lambda_{-})}{(t-1) (\lambda_{-}-\lambda_{+})} & \frac{t (\lambda_{-}-2 s)}{(t-1) (\lambda_{-}-\lambda_{+})} \\
 \frac{s t (\lambda_{+}-2 s t)}{(t-1) (\lambda_{-}-\lambda_{+})} & \frac{t (2 s-\lambda_{+})}{(t-1) (\lambda_{-}-\lambda_{+})} & \frac{s t (\lambda_{+}-2 s t)}{(t-1) (\lambda_{-}-\lambda_{+})} & \frac{t (2 s-\lambda_{+})}{(t-1) (\lambda_{-}-\lambda_{+})} \\
 \frac{1}{t-1} & \frac{1}{s-s t} & \frac{t}{t-1} & \frac{1}{s-s t} \\
 \frac{s t}{t-1} & \frac{t}{1-t} & \frac{s t}{t-1} & \frac{1}{1-t} \\
\end{array}
\right)
\end{align}
$$
