Is it possible to use numpy's linalg.matrix_power with a modulo so the elements don't grow larger than a certain value?
In order to prevent overflow, you can use the fact that you get the same result if you first take the modulo of each of your input numbers; in fact:
(M**k) mod p = ([M mod p]**k) mod p,
for a matrix M. This comes from the following two fundamental identities, which are valid for integers x and y (and a positive power p):
(x+y) mod p = ([x mod p]+[y mod p]) mod p # All additions can be done on numbers *modulo p*
(x*y) mod p = ([x mod p]*[y mod p]) mod p # All multiplications can be done on numbers *modulo p*
The same identities hold for matrices as well, since matrix addition and multiplication can be expressed through scalar addition and multiplication. With this, you only exponentiate small numbers (n mod p is generally much smaller than n) and are much less likely to get overflows. In NumPy, you would therefore simply do
((arr % p)**k) % p
in order to get (arr**k) mod p.
If this is still not enough (i.e., if there is a risk that [n mod p]**k causes overflow despite n mod p being small), you can break up the exponentiation into multiple exponentiations. The fundamental identities above yield
(n**[a+b]) mod p = ([{n mod p}**a mod p] * [{n mod p}**b mod p]) mod p
and
(n**[a*b]) mod p = ([n mod p]**a mod p)**b mod p.
Thus, you can break up power k as a+b+… or a*b*… or any combination thereof. The identities above allow you to perform only exponentiations of small numbers by small numbers, which greatly lowers the risk of integer overflows.