Matrix powers

Update:
 
GAUSS now has a function to compute matrix powers more efficiently. You can find the documentation here: powerM.
 

 
Original post:
 
The easiest way to calculate the power of a matrix is probably using the eigenvalues (which we will call lambda) and the eigenvectors (which we will call e_vec). The steps are:

1. Calculate the eigenvalues and eigenvectors:

   { lambda, e_vec } = eigv(A);

2. Diagonalize the eigenvalues:

   D = lambda .* eye(rows(A));

3. Multiply:

   power = 3;
   A_to_3 = e_vec * D.^power * inv(e_vec);

Here is a GAUSS procedure that with a simple example:

A = rndn(4, 4);
power = 3;
A_3 = A*A*A;
A_3_powmat = powMat(A, power);
print "A_3 = " A_3;
print "A_3_powmat = " A_3_powmat;

proc (1) = powMat(A, power);
    local e_val, e_vec, A_pow;
    
    //Calculate eigenvalues
    { e_val, e_vec } = eigv(A);
    
    //Calculate power of matrix
    A_pow = e_vec * (e_val .* eye(rows(A))).^power * inv(e_vec);
    
    //The eigenvalues for many matrices are complex
    //and will return a complex portion that is round-off error
    //remove the round-off error for real matrices
    if iscplx(A);
        retp(A_pow);
    else;
        retp(real(A_pow));
    endif;
endp;

Thanks a lot for the answer, that is one possible way to calculate it. After digging into the language reference, I have found a command that seems even easier.
For example, if I want to calculate the nth power of a matrix A, I can do the following:

A_pow = polyeval(A, 1|zeros(n,1))

I do not know how this function is evaluated, so I can't say anything about how round-off errors are accumulated with the two methods. If the matrix contains very small numbers, the eigenvalue decomposition and polyeval can lead to different results, but no idea which is more reliable.

Thanks again for your help.

Szabolcs