Exact powers of an integer matrix Assume that I have a large, non-symmetric matrix $\bf{A}$ of zeros and ones. If I want the exact diagonal elements of $\bf{A}^{k}$ (not the trace, but the elements themselves) where $k$ is some large whole number, I can


*

*Multiply the matrices out using with exact integers, possibly speeding the computation by using some Addition-chain exponentiation.

*Diagonalize the matrix as $\bf A = V \Gamma V^{-1}$ and compute $\bf A^{k} = V \Gamma^{k} V^{-1}$
The first method is slow but accurate, while the second method can introduce numerical errors (since the exact answer is required). Is there some decomposition I am missing that does what I'm looking for?
 A: (edit : this is for large k, so it may not be useful for k comparable to n)
A faster way to compute high powers of a matrix is by computing in the $\mathbb{Z}$-algebra $\mathbb{Z}[M]$.
The space $(I_n,M,M^2, \ldots)$ is a subspace of $M_n(\mathbb{Z})$, and in fact the Cayley-Hamilton theorem tells you that it has dimension at most $n$. For example, sometimes, the $n$ coefficients on the diagonal of the matrix (or the the coefficients on the first row, or column, etc) are enough to know what the rest of the coefficients are.
If $P$ is the minimal polynomial of $M$ over $\mathbb{Z}$ and has degree $d = d_P$, the basis the best suited in order to compute products, however, is the basis $(1,M,\ldots M^{d-1})$ : we want to use the isomorphism $\mathbb{Z}[M] = \mathbb{Z}[X]/P$.
Doing one multiplication in there requires $O(d^2)$ operations :
Keep a table of the coordinates of $(1,M,\ldots M^{2d-2})$ in the basis $(1,M,\ldots M^{d-1})$, then compute the coefficients $c_k$ in
$$ (\Sigma_{i=0}^{d-1} a_i M^i)(\Sigma_{j=0}^{d-1} b_j M^j) = \Sigma_{k=0}^{2d-2} (c_k = \Sigma_{i+j=k} a_i b_j) M^k$$
Then use that table to convert this sum into a vector $\Sigma_{k=0}^{d-1} d_k M^k$.
You can further optimize this algorithm if the polynomial $P$ is reducible.
If $P = QR$ with $Q$ and $R$ having no factor in common, then $\mathbb{Z}[M] = \mathbb{Z}[X]/P \simeq \mathbb{Z}[X]/Q \times \mathbb{Z}[X]/R$. The isomorphism here sends $X$ onto the pair $X,X$. you will need to compute the inverse of the isomorphism, but once you do that, you gain in complexity because $d_Q^2 + d_R^2 < (d_Q + d_R)^2 = d_P^2$.
You can also improve if you have $P = Q^k$ for some $k$, because then, $\mathbb{Z}[X]/P \simeq \mathbb{Z}[X,Y]/(Q(X),Y^k)$, where the isomorphism sends $X$ onto $X+Y$. Once again, the computations there are shorter since you have $k$ polynomials of degree $d/k$ to compute, it will take $O(k(d/k)^2 = O(d^2/k)$.
You have to compute $P$, its factorization, and the isomorphisms only once. Once you have done all this preliminary work you get a very good algorithm. 
If the matrix is diagonalizable (in $M_n(\mathbb{Z})$), you will get that $P$ factors completely into distincts factors of degree 1, and this will give you an algorithm equivalent to the one you gave.
Even if you don't want to factorize $P$, simply working in $\mathbb{Z}[X]/P$ gives a good optimisation. You can compute $P$ by taking some element $x$ in $\mathbb{Z}^n$, and checking when $M^k x$ becomes linearly dependant on $1, x, \ldots x^{k-1}$. If you are lucky you get some dependancy for $k<n$, which means you have found a factor of $P$ of degree $k<n$. If you know of any subspace of $\mathbb{Z}^n$ stable by $M$, be sure to make use of it, picking $x$ in it will give you such a factor of $P$.
