@Carl Love 

Problem exactly happened when I ran program Matrix_powers_finite_field.mv (at http://www.mapleprimes.com/view.aspx?sf=215285_Answer/Matrix_powers_finite_field.mw)with modulus polynomial: Z^8+Z^7+Z^6+Z+1 and Matrix:

M:<x^7+x^6+x^5+x^2+1, x^7+x^6+x^4+x^3+x^2+x, x^6+x^5+x^4+x^2, x^7+x^6+x^3+x+1;

x^3+x^5, x^6+x^5+x^3+1, x^7+1, x^6+x^4+x^2+1;

x^6+x^2+x, x^7+x^6+x^5+x^3+x^2, x^5+x+1, x^6+x^5+x^4+x^3+x^2+1;

x^6, x^6+x^5+x^4+x^3+x, x^6+x^5+x^4+x^2, x^7+x^4+x >;

That is: 

- Extract the nonlinear irreducible parts (this command works even if there is more than one irreducible factor):

> irr:= map(collect, select(q-> degree(q,t) > 1, map2(op,1,CF[2])), t);

> encode(irr,x) --> t^3 +209 t^2 + 11 t + 223;

>q:= irr[1]:

Now construct the splitting field that contains q's roots.

>alias(z= RootOf(q)):

- Command ">PF:= Primfield({x,z}) mod p:" don't run

Perhaps the degree of polynomial irr (and q) too large

Modulus polynomial: Z^8+Z^7+Z^6+Z+1

Matrix: 

M:<x^7+x^6+x^5+x^2+1, x^7+x^6+x^4+x^3+x^2+x, x^6+x^5+x^4+x^2, x^7+x^6+x^3+x+1;

x^3+x^5, x^6+x^5+x^3+1, x^7+1, x^6+x^4+x^2+1;

x^6+x^2+x, x^7+x^6+x^5+x^3+x^2, x^5+x+1, x^6+x^5+x^4+x^3+x^2+1;

x^6, x^6+x^5+x^4+x^3+x, x^6+x^5+x^4+x^2, x^7+x^4+x >;

OR

Modulus polynomial: Z^8+Z^5+Z^3+Z+1

Matrix: 

M:<x^3+x^4+x^6+x^7+x, 1+x^2+x^4+x^5+x^6, 1+x+x^2+x^3, x^7+x^6+x^5+x^4;

x^7+x^5+x^4+x^3, x^6+x^4+x^2+1, x^3+x^4+x^6+x^7+x^2+x+1, 1+x^2+x^3+x^4+x^5;

x^7+x^5+x^2, x^7+x^5+x^3+x^2+1, x+x^2+x^6, x^2+x^3+x^5;

x^3+x^4+x^6+1, 1+x^2+x^3+x^4, x^6+x^5+x^4+x^3, x^7+x^3 >;

Thanks alot

@Carl Love Thank you very much

@Carl Love 

I want to compute rank of matrix M + I . But no result. This is used command

> A:=Expand~(M+I) mod p; {I: IdentityMatrix}

> Gausselim( A ,'rank') mod p ;

Error, (in mod/Gausselim) 1st argument must be a matrix over a finite field

Please help me! Thank you very much.

@Carl Love 

- In the matrix M, M[i,j] ≠ 0 (0<i,j<=n)

- For item (3) not divide V[i] ( divide Sum², i.e M[i,j]/Sum²)

@Carl Love 

- The finite field and matrix in this question are the same as in question "How to find k^th root of the given matrix over finite field (at URL: http://www.mapleprimes.com/questions/203997-How-To-Find-Kth-Root-Of-The-Given-Matrix#comment215683)".

- Dividing by row sum squared M[i,j]/V[i] <=> M[i,j] x (V[i])^-1 ((V[i])^-1 is inverse of V[i])

@Carl Love 

Could You help me again to compute direct exponentiation of the given matrix M in finite field?

(If M = [m_i,j], we say Direct Exponent (element-wise exponent matrix), M_d^k of M is a matrix whose each element is the result of exponentiation of corresponding elements of M. If k=2, then we say M_d² is a direct square matrix of M )

@Carl Love

I'm sorry!! I have mistakenly provided the matrixs M, R.

True Matrixs:  

R= [[0,1,0,0],[0,0,1,0],[0,0,0,1],[x^2,1,x,x]

M=[[4,1,2,2],[8,6,5,6],[11,14,10,9],[2,2,15,11]]

I successfully implemented Your code. Thank You very much

@Carl Love 

I'm sorry!! I have mistakenly provided the matrix R. Matrix true R= [[0,1,0,0],[0,0,1,0],[0,0,0,1],[x^2,x,1,1]. With this R, R^4 = M.

@Carl Love 

Of course I known R^4 = M and I checked by simply comnand Expand~(R^4) mod p;. But the problem is that if I excuted Your proc with M:= [ [x^2, x,1,1], [x^3,x^2+x,x^2+1,x^2+x],[x^3+x+1,x^3+x^2+x, x^3+x,x^2+x],[x,x,x^3+x^2+x+1, x^3+x+1], but not receive R= [[0,1,0,0],[0,0,1,0],[0,0,0,1],[1,1,x,x^2].

And with the following matrix

M:= < x^3+x^2+x, x+1, x^3+x+1, x^5+x+1; 
 x, x, 1, x+1; 
 x^3+x^2+1, 1, x^3+x+1, x^5+x^3; 
 x^4+x^3+x, x^2+x, x^4+x^2+1, x^6+x^2+x+1 >:

also not receive R=[[0,1,0,x],[1,0,1,0],[x,0,1,3],[1,0,1,x^2]], while R^3 = M {Expand~(R^3) mod p = M in GF(2^8)/(z^8+z^4+z^3+z+1)}.

Moreover, the input of this proc is any matrix, need to find R (R is completely unknown)

@Carl Love

Thank You very much. I successfully implemented your code with the given matrix above. But with different matrix

M:= < x^2,           x,               1,                   1;
         x^3,           x^2+x,        x^2+1,            x^2+x;
         x^3+x+1,    x^3+x^2+x, x^3+x,            x^2+x;
         x,               x,               x^3+x^2+x+1, x^3+x+1 >:

(or this matrix: 

M:= < x^3+x^2+x, x+1, x^3+x+1, x^5+x+1;
 x, x, 1, x+1;
 x^3+x^2+1, 1, x^3+x+1, x^5+x^3;
 x^4+x^3+x, x^2+x, x^4+x^2+1, x^6+x^2+x+1 >:

)

There is error with this matrix when applying the isomorphism to R

> R:= Iso~(R);

Error, invalid left hand side in assignment

- In this case R will be the following matrix:

R = <0, 1, 0, 0 ;

        0, 0, 1, 0;

        0, 0, 0, 1:

        1, 1, x, x^2>

(In GF(2^8)/(z^8 + z^4 + z^3 +z +1) R^4 = M)

Could You help me again to check with this matrix M ?  

@Carl Love 

Please help me convert the entries of R (and R^4) from GF(2^16) back to GF(2^8). 

{ I think, GF(2^8) is the sub set of GF(2^16). So Diag[1,1] and Diag[2,2] of the new matrix newM are not changed, ie in matrix newM Diag[1,1] =37 and Diag[2,2] =217). But in our results, these values have been changed completely. And while convert R and R^4 back to GF(2^8), the entries ≤ 255 will be not changed}

Furthermore, for the calculation of the roots of the factor (t² +231t +30), can You choose an irreducible polynomial ext2 over GF(2^8)? (ie ext2 = t² + at + b, where a, b in GF(2^8). GF((2⁸)² is the composition field)

@Carl Love

In your answer, R⁴ = newM  and some elements of R in GF(2^16).

@Carl Love 

- It's right that, I need to explicitly raise a matrix (Diag)to the 1/4 power (as sqrt(sqrt(Diag)) ).

- In my project, Diag[1,1] = 37 (this is root of factor t + 37), Diag[2,2]=217 (this is root of t + 217)and Diag[3,3] = x₃, Diag[4,4]=x₄ , where x3 and x4 are roots of t² + 231t +30, these roots in GF(2⁸)² 

- Need to compute R = P x Diag^^(1/4) X P^-1 . (Note 1/4 power, not 4 power)

(Example: for matrix M above, R = [[0,1,0,0],[0,0,1,0],[0,0,0,1],[1,2,1,4]] (Companion matrix) and in "GF(2^8)/...." R^4 = M).

@Carl Love 

Thank you very much

Your answer nearly meets my needs, it executes on my version of Maple. Exists only afew small problems as follows:

1) The matrix D is not take powers 1/4. (D = Diag^1/4 ).

2) Then the matrix R = P x D x P^-1 . Then elements of the received R will return to GF(2⁸) and R⁴ = M (The calculation of R⁴ is performed in GF(2⁸)/(y^8 + y^4 + Y^3 + y +1))