Computing the Frobenius form of a matrix

[Bill Allombert <allomber@math.u-bordeaux.fr>]

Hello PARI developers,

The attached script implement matrix reductions algorithms:

matfrobenius(M{,flags}): compute the Frobenius form F of matrix M. If
flag is 1, return [F,P] where P is the basis change.

matbasischange(M,N}): Compute the basis change matrix P so that
N=P*M*P^-1. Return 0 if such matrix does not exist.
 
The algorithm use SNF computation over k[X] so work better with exact
coefficients, but is generic.

Using matbasischange() and matfrobenius(), it is easy to implement
Jordan reduction even using POLMOD for the eigen value instead of
complex approximation.

I could add a function for computing elementary divisors of a 
matrix but it is trivial: 
matelemdivisors(M)=matsnf(M-'x*matid(#M),6)

Before adding thoses features to PARI, I would like to hear comments
about the algorithms and the terminology (how to name the functions,
etc.).

Cheers,
Bill.

build_basischange(N,A,B)=
{
        local(n,U);
        n=#N;
        U=B[1]^-1*A[1];
        concat(vector(n,j,
               Mat(sum(i=1,n,
                        subst(U[i,j],'x,N)[,i])
                        )));
}

matbasischange(M,N)=
{
        local(A,B,U,R,n);
        n=#M;
        A=matsnf(M-'x*matid(n),3);
        B=matsnf(N-'x*matid(n),3);
        if(A[3]!=B[3],print("warning: not semblable:",A[3]-B[3]));
        build_basischange(N,A,B);
}

addhelp(matbasischange,"matbasischange(M,N}): Compute the basis change matrix P so that N=P*M*P^-1. Return 0 if such matrix does not exist.");

build_matfrobenius(n,R)=
{ 
  local(P,d,s);
  concat(vector(#R,i,
    P=R[i];d=poldegree(P);s+=d;
    concat([matrix(d,s-d),matcompanion(P),matrix(d,n-s)])~
  ))~  
}

matfrobenius(M,flags)=
{
  local(n,A,B,N);
  n=#M;
  if (flags==0,
    A=matsnf(M-'x*matid(n),6);
    build_matfrobenius(n,A);
  ,
    A=matsnf(M-'x*matid(n),3);
    N=build_matfrobenius(n,A[3]);
    B=matsnf(N-'x*matid(n),3);
    R=build_basischange(N,A,B);
    [N,R])
}
addhelp(matfrobenius,"matfrobenius(M{,flags}): compute the Frobenius form F of matrix M. If flag is 1, return [F,P] where P is the basis change.");