Bill Allombert on Sun, 13 Jun 2004 19:30:26 +0200

 Computing the Frobenius form of a matrix

```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.");
```