share/matrix/rateig.mac

   1 /*-*-macsyma-*-*/
   2 base(pol):=if atom(pol) or inpart(pol,0) # "^" then pol else inpart(pol,1)$
   3
   4 rateigenvectors(mat,[var]):=
   5      block([scalarmatrixp:false,
   6             listarith:algebraic:inflag:true,
   7             ans:[],
   8             tellrats:tellrat(),
   9             id:ident(length(mat))],
  10            var:if var=[] then 'lambda else var[1],
  11            for pol in 'false*factor(charpoly(mat,var)) do
  12             if not(freeof(var,pol))
  13              then (rat(base(pol),var),
  14                    tellrat(%%/coeff(%%,var,hipow(%%,var))),
  15                    ans:cons(cons(pol,args(transpose(%rateig(echelon(mat-var*id),id)))),
  16                             ans)),
  17            apply('tellrat,tellrats),
  18            ans)$
  19
  20 %rateig(mat,id):=
  21       if mat=-mat then id
  22  else if mat*id=id then zeromatrix(length(mat),0)
  23  else if mat[1,1]=0 then addcol(col(id,1),
  24                                  (id:submatrix(id,1))
  25                                 .
  26                                  %rateig(submatrix(rest(mat,-1),1),
  27                                          rest(id)))
  28  else  addrow(matrix(-rest(mat[1])),id:submatrix(1,id,1))
  29       .
  30        %rateig(submatrix(1,mat,1),id)$