murraymiller:=proc(p,m)
  local l:=Dimension(p)[1],P,s,i,z,t,T,r,c:
  P:=p:
  for s from 1 to l do:
    z:=[seq(P[s,i],i=(s+1)..l)]:
    if andmap(temp->temp=0, z) then: return [s,P]: fi:
    if P[s,s+1]=0 then:
      for t from s+2 to l do:
        if P[s,t]<>0 then: break: fi:
      od:
      P:=RowOperation(P,[s+1,t]):
      P:=ColumnOperation(P,[s+1,t]):
    fi: 
    T:=Matrix(l,l,(r,c)-> `if`(r=c,1,0)):
    T[[s+1],(s+1)..l]:=Matrix([ [seq(P[s,i],i=(s+1..l))] ]):
    P:=simplify(subs(x=x*q^(-m),T).P.MatrixInverse(T)):
  od:  
end:

mmrec:=proc(p,m,var)
  local sf, sb, i, v, f, temp, gf, gb, fin:
  local l:=Dimension(p)[1], gseq:=[seq(g(x/q^(m*t)),t=-l-2..l+2)]:
  gf:=fun->subs(x=x*q^m,fun):
  gb:=fun->subs(x=x*q^(-m),fun):

  for i from 1 to l-1 do:
    v:=convert( [seq( sf(f[i]),i=1..l)],Vector):
    temp:=-(p[i].v - sf(f[i+1]))+f[i]:
    f[i+1]:=sb(temp):
    f[i+1]:=subs(f[1]=g(x),f[i+1]);
    f[i+1]:=expand(eval(subs(sf=gf,sb=gb,f[i+1])));
    f[i+1]:=collect(f[i+1],gseq,simplify);
  od:
  
  v:=convert( [seq( sf(f[i]),i=1..l)],Vector):
  fin:= p[l].v - f[l]:

  fin:=subs(f[1]=var(x),fin):

  fin:=expand(eval(subs(sf=gf,sb=gb,fin))):
  fin:=collect(fin,gseq,simplify):
  
  return(fin):
end:

tt:=Matrix(7, 7, [[1,x*q^2,x^2*q^4,x*q,x^2*q^2,0,0],[0,x*q^2,0,0,0,1,0],[0,0,0,0,0,0,1],[0,x*q^2,0,x*q,0,1,0],[0,0,0,0,x*q^2,0,1],[1,x*q^2,x^2*q^4,x*q,0,0,0],[1,x*q^2,x^2*q^4,0,0,0,0]]):