Compute matrix power with integer entries

Question

Let's say I have a constant matrix A and I want to compute pow(A, n). As described in this question I can calculate its eigenvalue decomposition (or more generally, its invariant subspaces and the generalized modal matrix) to speed up the process.

If A is a square matrix of size k, then the algorithm has complexity O(k log n) via exponentiation by squaring, and a preparation cost (to compute the modal matrix) of O(k^3).

The problem I am thinking about is loss of precision. Calculating eigenvalues et al takes us out of the domain of integers into floating point numbers. Even though in the end, we know that pow(A, n) has to have all integer entries, the algorithm outlined above only computes floating point numbers.

Another way is to exploit only exponentiation by squaring but that gives us only a O(k^3 log n) algorithm.

Is there a way to accurately - without converting to floating point numbers - compute pow(A, n) fast?

Answer 1

Eigenvalue decomposition is also possible for a matrix over a finite field, but only if the field is just right. So it not just takes preprocessing to do the eigenvalue decomposition, but also to find (some) finite field(s) over which that is even possible.

Finding multiple solutions is useful to avoid having work with gigantic finite fields, then compute pow(A, n) in some small fields and use the CRT to work out what the solution would have been in ℤ. But this requires somehow having a sufficient number of fields of sufficient size to work with and you wouldn't really know in advance what will be sufficient (there is always some n above which it stops working), so maybe this all won't work in practice.

As a small example, take:

A = [[1, 1],
     [1, 0]]

Characteristic x² - x - 1, let's guess that modulo 1009 will work (it does), then there are roots 383 and 627, so:

A = QDP mod 1009
Q = [[  1,   1],
     [382, 626]]
D = [[383,   0],
     [  0, 627]]
P = [[ 77, 153],
     [933, 856]]

So for example

pow(A, 15) = Q [[928,   0], P = [[987, 610],
                [  0, 436]]      [610, 377]]

Fibonacci numbers as expected, so it all worked out. But with just 1009 as the modulus, going above 15 for the exponent makes the result not match would it would be in ℤ, then we would need more/bigger fields.