One missing maxima feature, that I have the math behind, but not a program

Updated the project at https://github.com/Mickle-Mouse/Maxima-Matrix-Power


Below is also a copy of the code and a few simple test cases:
/* Michael Valenzuela's original symbolic matrix power */
/* This very non-robust code is effectively released   */
/* Under the "New BSD License"                         */
/* This is compatable with GPL, but is not copy-left   */
/* I just don't want to be held accountable for this   */
/* code                                                */


/* Version 1.1 --- non robust, not very effecient version */
matrix_power(A, t):=block([A_, A_list, B_, B_list, eqs, eigvals, I_, i, j,
k, numeq, tt, solvedeq, offset],
/* A should be a matrix          */
/* t should be a scalar variable */

/* Initialize Variables */
declare(A_, nonscalar),
declare(I_, nonscalar),
/* This is the value we set t to this value. Should be >=0 */
offset: 0,
/* The number of equations */
numeq: length(A),
/* The list of relevant variables */
eqs: makelist(0,i,1,numeq),
A_list: makelist( A_^^(i+offset)=tt,i,0,numeq-1),
B_list: makelist( B_[i],i,1,numeq),


/* Correct for A_^^0 being 1, it should be identity */
if (offset=0) then
  A_list[1]: I_=tt else
  0,


/* Exact Eigenvalues and calculate powers of A */
eigvals:rat(radcan(fullratsimp(eivals(A)))),
for i:1 thru numeq do block([],
  A_list[i]: at( A_list[i], tt=A^^(i+offset-1))
),
A_list : reverse(A_list),


/* Setup the first equation */
k:1,
eqs[1]: 0,
for i:1 thru length(eigvals[1]) do block([],
  for j:1 thru eigvals[2][i] do block ([],
    if( eigvals[1][i] # 0 ) then
      eqs[1]: eqs[1] + B_list[k] * t^(j-1) * eigvals[1][i]^t else
      eqs[1]: eqs[1] + B_list[k] * kron_delta(j-1,t),
    k:k+1
  )
),
eqs[1]: A_^^t = eqs[1],


/* Generate all the necessary equations */
for i:2 thru numeq do block([],
  eqs[i]: at(eqs[i-1], [t=1+tt]),
  eqs[i]: at(eqs[i], [tt=t])
),


/* Lets take care of the case where A_^^0 turns into 1 */
if (offset = 0) then
  eqs[1]: at(eqs[1], A_^^t=I_) else
  0,


/* Set t=offset and solve the equations */
solvedeq: linsolve( at(eqs, t=offset), B_list ) ,


/* Substitute the symbolic A_^k with the calculation */
solvedeq: radcan( at(solvedeq, A_list)),


/* Substitute the resulting equations into eq0 */
return( rhs(at(eqs[1], solvedeq ) ) )
);

m1: matrix([1, 1, 0, 0],[0, 1, 1, 0],[0,0,1,-1/8],[0,0,1/2,1/2])$
m2: matrix([1, 0, 1],[0, 1, 0],[0,0,0])$
shift_matrix: matrix([0, 1, 0, 0],[0, 0, 1, 0],[0,0,0,1],[0,0,0,0])$

m1power(t):= (''(matrix_power(m1, t)));
m2power(t):= (''(matrix_power(m2, t)));
shift_matrixpower(t):= (''(matrix_power(shift_matrix, t)));

for i:-10 thru 20 do
    display(is( m1^^i = m1power(i)));

for i:0 thru 20 do
    display(is( m2^^i = m2power(i)));

is( shift_matrix^^0 = shift_matrixpower(0));
is( shift_matrix^^1 = shift_matrixpower(1));
is( shift_matrix^^2 = shift_matrixpower(2));
is( shift_matrix^^3 = shift_matrixpower(3));
is( shift_matrix^^4 = shift_matrixpower(4));




On Sat, Dec 21, 2013 at 4:32 PM, Mike Valenzuela <mickle.mouse at gmail.com>wrote:

> Here is a working example:
>
> m1: matrix([1, 1, 0, 0],[0, 1, 1, 0],[0,0,1,-1/8],[0,0,1/2,1/2])$
> sympower(t):= (''(matrix_power(m1, t)));
>
> %o
>
> sympower(t):=matrix([1,t,4^(2-t)*3^(t+1)+t*4^(2-t)*3^(t-1)+8*t-48,-4^(2-t)*3^t-2*t*4^(1-t)*3^(t-1)-2*t+16],[0,1,-2*4^(1-t)*3^t-t*4^(1-t)*3^(t-1)+8,(2*3^t)/4^t+(2*t*3^(t-1))/4^t-2],[0,0,3^t/4^t+(t*3^(t-1))/4^t,-(t*3^(t-1))/(2*4^t)],[0,0,(2*t*3^(t-1))/4^t,3^t/4^t-(t*3^(t-1))/4^t])
>
>
> is( m1^^15 = sympower(15));
> is( m1^^-1 = sympower(-1));
>
> %o
> True
> %o
> True
>
>
>
>
>
> On Sat, Dec 21, 2013 at 8:21 AM, Barton Willis <willisb at unk.edu> wrote:
>
>>  Two more examples:
>>
>> (%i14) m : matrix([5,7,1],[9,11,7],[5,7,8])$
>>
>> OK, but suboptimal (linsolve detects under determined equations)
>>
>> (%i23) matrix_power(m,1);
>> solve: dependent equations eliminated: (3)
>> (%o23)
>> matrix([-5*(4*%r5-1)+(sqrt(2)*sqrt(3)*sqrt(19)+10)*%r5+(10-sqrt(2)*sqrt(3)*sqrt(19))*%r5,(sqrt(2)*sqrt(3)*sqrt(19)+10)*%r5+(10-sqrt(2)*sqrt(3)*sqrt(19))*%r5-20*%r5+7,(sqrt(2)*sqrt(3)*sqrt(19)+10)*%r5+(10-sqrt(2)*sqrt(3)*sqrt(19))*%r5-20*%r5+1],[(sqrt(2)*sqrt(3)*sqrt(19)+10)*%r5+(10-sqrt(2)*sqrt(3)*sqrt(19))*%r5-20*%r5+9,(sqrt(2)*sqrt(3)*sqrt(19)+10)*%r5+(10-sqrt(2)*sqrt(3)*sqrt(19))*%r5-20*%r5+11,(sqrt(2)*sqrt(3)*sqrt(19)+10)*%r5+(10-sqrt(2)*sqrt(3)*sqrt(19))*%r5-20*%r5+7],[-5*(4*%r5-1)+(sqrt(2)*sqrt(3)*sqrt(19)+10)*%r5+(10-sqrt(2)*sqrt(3)*sqrt(19))*%r5,(sqrt(2)*sqrt(3)*sqrt(19)+10)*%r5+(10-sqrt(2)*sqrt(3)*sqrt(19))*%r5-20*%r5+7,(sqrt(2)*sqrt(3)*sqrt(19)+10)*%r5+(10-sqrt(2)*sqrt(3)*sqrt(19))*%r5+4*(2-5*%r5)])
>>
>> (%i24) ratsimp(%);
>> (%o24) matrix([5,7,1],[9,11,7],[5,7,8])
>>
>> Likely completely bogus:
>>
>> (%i30) ratsimp(matrix_power(m,1/2). matrix_power(m,1/2));
>> solve: dependent equations eliminated: (3)
>> solve: dependent equations eliminated: (3)
>>  << Expression too long to display! >>
>>
>> OK:
>>
>> (%i31) ratsimp(matrixfun(lambda([x], sqrt(x)), m) . matrixfun(lambda([x],
>> sqrt(x)), m));
>> (%o31) matrix([5,7,1],[9,11,7],[5,7,8])
>>
>>
>>
>>  --Barton
>>
>>
>>  ------------------------------
>> *From:* maxima-bounces at math.utexas.edu <maxima-bounces at math.utexas.edu>
>> on behalf of Barton Willis <willisb at unk.edu>
>> *Sent:* Saturday, December 21, 2013 07:21
>> *To:* Mike Valenzuela
>> *Cc:* maxima at math.utexas.edu
>> *Subject:* Re: [Maxima] One missing maxima feature, that I have the math
>> behind, but not a program
>>
>>
>> Examples:
>>
>>
>>  (%i2) m : matrix([5,7],[9,11])$
>>
>>
>>  Suboptimal?
>>
>>
>>      (%i3) matrix_power(m,0);
>>       expt: undefined: 0^0
>>
>>       0: matrix_power(a=matrix([5,7],[9,11]),t=0)(mp.mac line 35)
>>
>>
>>  For the 0-th power, matrixfun gives the identity:
>>
>>
>>    (%i4) matrixfun(lambda([x], x^0), m);
>>    (%o4) matrix([1,0],[0,1])
>>
>>
>>  For an inverse, matrix_power returns a scalar
>>
>>
>>     (%i5) matrix_power(m,-1);
>>    (%o5) B_[2]/(3*2^(3/2)+8)+B_[1]/(8-3*2^(3/2))
>>
>>
>>     (%i6) matrixfun(lambda([x], x^-1), m);
>>     (%o6) matrix([-11/8,7/8],[9/8,-5/8])
>>
>>
>>      (%i9) m : matrix([a,b],[b,c]);
>>     (%o9) matrix([a,b],[b,c])
>>
>>
>>  Symbolic matrices:
>>
>>
>>      (%i18) m : matrix([a,b],[b,c]);
>>     (%o18) matrix([a,b],[b,c])
>>
>>
>>      (%i19) matrix_power(m,2);
>>     (%o19)
>> (B_[2]*(sqrt(c^2-2*a*c+4*b^2+a^2)+c+a)^2)/4+(B_[1]*(sqrt(c^2-2*a*c+4*b^2+a^2)-c-a)^2)/4
>>
>>
>>
>>     (%i20) matrixfun(lambda([x],x^2),m);
>>     Proviso: assuming c^2-2*a*c+4*b^2+a^2 # 0
>>     (%o20) matrix([b^2+a^2,b*c+a*b],[b*c+a*b,c^2+b^2])
>>
>>
>>   --Barton (author of matrixfun)
>>
>>
>>  ------------------------------
>>
>
>