//////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Linear Algebra";
info="
LIBRARY:  linalg.lib  Algorithmic Linear Algebra
AUTHORS:  Ivor Saynisch (ivs@math.tu-cottbus.de)
@*        Mathias Schulze (mschulze@mathematik.uni-kl.de)

PROCEDURES:
 inverse(A);          matrix, the inverse of A
 inverse_B(A);        list(matrix Inv,poly p),Inv*A=p*En ( using busadj(A) )
 inverse_L(A);        list(matrix Inv,poly p),Inv*A=p*En ( using lift )
 sym_gauss(A);        symmetric gaussian algorithm
 orthogonalize(A);    Gram-Schmidt orthogonalization
 diag_test(A);        test whether A can be diagnolized
 busadj(A);           coefficients of Adj(E*t-A) and coefficients of det(E*t-A)
 charpoly(A,v);       characteristic polynomial of A ( using busadj(A) )
 adjoint(A);          adjoint of A ( using busadj(A) )
 det_B(A);            determinant of A ( using busadj(A) )
 gaussred(A);         gaussian reduction: P*A=U*S, S a row reduced form of A
 gaussred_pivot(A);   gaussian reduction: P*A=U*S, uses row pivoting
 gauss_nf(A);         gaussian normal form of A
 mat_rk(A);           rank of constant matrix A
 U_D_O(A);            P*A=U*D*O, P,D,U,O=permutaion,diag,lower-,upper-triang
 pos_def(A,i);        test symmetric matrix for positive definiteness
 hessenberg(M);       Hessenberg form of M
 eigenvals(M);        eigenvalues with multiplicities of M
 minipoly(M);         minimal polynomial of M
 spnf(sp);            normal form of spectrum sp
 spprint(sp);         print spectrum sp
 jordan(M);           Jordan data of M
 jordanbasis(M);      Jordan basis and weight filtration of M
 jordanmatrix(jd);    Jordan matrix with Jordan data jd
 jordannf(M);         Jordan normal form of M
";

LIB "matrix.lib";
LIB "ring.lib";
LIB "elim.lib";
LIB "general.lib";
//////////////////////////////////////////////////////////////////////////////
// help functions
//////////////////////////////////////////////////////////////////////////////
static proc const_mat(matrix A)
"RETURN:   1 (0) if A is (is not) a constant matrix"
{
  int i;
  int n=ncols(A);
  def BR=basering;
  def @R=changeord("dp,c",BR);
  setring @R;
  matrix A=fetch(BR,A);
  for(i=1;i<=n;i=i+1){
    if(deg(lead(A)[i])>=1){
      //"input is not a constant matrix";
      kill @R;
      setring BR;
      return(0);
    }
  }
  kill @R;
  setring BR;
  return(1);
}
//////////////////////////////////////////////////////////////////////////////
static proc red(matrix A,int i,int j)
"USAGE:    red(A,i,j);  A = constant matrix
          reduces column j with respect to A[i,i] and column i
          reduces row j with respect to A[i,i] and row i
RETURN:   matrix
"
{
  module m=module(A);

  if(A[i,i]==0){
    m[i]=m[i]+m[j];
    m=module(transpose(matrix(m)));
    m[i]=m[i]+m[j];
    m=module(transpose(matrix(m)));
  }

  A=matrix(m);
  m[j]=m[j]-(A[i,j]/A[i,i])*m[i];
  m=module(transpose(matrix(m)));
  m[j]=m[j]-(A[i,j]/A[i,i])*m[i];
  m=module(transpose(matrix(m)));

  return(matrix(m));
}

//////////////////////////////////////////////////////////////////////////////
proc inner_product(vector v1,vector v2)
"RETURN:   inner product <v1,v2> "
{
  int k;
  if (nrows(v2)>nrows(v1)) { k=nrows(v2); } else { k=nrows(v1); }
  return ((transpose(matrix(v1,k,1))*matrix(v2,k,1))[1,1]);
}

/////////////////////////////////////////////////////////////////////////////
// user functions
/////////////////////////////////////////////////////////////////////////////

proc inverse(matrix A, list #)
"USAGE:    inverse(A [,opt]);  A a square matrix, opt integer
RETURN:
@format
          a matrix:
          - the inverse matrix of A, if A is invertible;
          - the 1x1 0-matrix if A is not invertible (in the polynomial ring!).
          There are the following options:
          - opt=0 or not given: heuristically best option from below
          - opt=1 : apply std to (transpose(E,A)), ordering (C,dp).
          - opt=2 : apply interred (transpose(E,A)), ordering (C,dp).
          - opt=3 : apply lift(A,E), ordering (C,dp).
@end format
NOTE:     parameters and minpoly are allowed; opt=2 is only correct for
          matrices with entries in a field
SEE ALSO: inverse_B, inverse_L
EXAMPLE:  example inverse; shows an example
"
{
//--------------------------- initialization and check ------------------------
   int ii,u,notInvertible,opt;
   matrix invA;
   int db = printlevel-voice+3;      //used for comments
   def R=basering;
   string mp = string(minpoly);
   int n = nrows(A);
   module M = A;
   intvec v = option(get);           //get options to reset it later
   if ( ncols(A)!=n )
   {
     ERROR("// ** no square matrix");
   }
//----------------------- choose heurisitically best option ------------------
// This may change later, depending on improvements of the implemantation
// at the monent we use if opt=0 or opt not given:
// opt = 1 (std) for everything
// opt = 2 (interred) for nothing, NOTE: interred is ok for constant matricea
// opt = 3 (lift) for nothing
// NOTE: interred is ok for constant matrices only (Beispiele am Ende der lib)
   if(size(#) != 0) {opt = #[1];}
   if(opt == 0)
   {
      if(npars(R) == 0)                       //no parameters
      {
         if( size(ideal(A-jet(A,0))) == 0 )   //constant matrix
         {opt = 1;}
         else
         {opt = 1;}
      }
      else {opt = 1;}
   }
//------------------------- change ring if necessary -------------------------
   if( ordstr(R) != "C,dp(nvars(R))" )
   {
     u=1;
     def @R=changeord("C,dp",R);
     setring @R;
     module M = fetch(R,M);
     execute("minpoly="+mp+";");
   }
//----------------------------- opt=3: use lift ------------------------------
   if( opt==3 )
   {
      module D2;
      D2 = lift(M,freemodule(n));
      if (size(ideal(D2))==0)
      {                                               //catch error in lift
         dbprint(db,"// ** matrix is not invertible");
         setring R;
         if (u==1) { kill @R;}
         return(invA);
      }
   }
//-------------- opt = 1 resp. opt = 2: use std resp. interred --------------
   if( opt==1 or opt==2 )
   {
      option(redSB);
      module B = freemodule(n),M;
      if(opt == 2)
      {
         module D = interred(transpose(B));
         D = transpose(simplify(D,1));
      }
      if(opt == 1)
      {
         module D = std(transpose(B));
         D = transpose(simplify(D,1));
      }
      module D2 = D[1..n];
      module D1 = D[n+1..2*n];
//----------------------- check if matrix is invertible ----------------------
      for (ii=1; ii<=n; ii++)
      {
         if ( D1[ii] != gen(ii) )
         {
            notInvertible = 1;
            break;
         }
      }
   }
   option(set,v);
//------------------ return to basering and return result ---------------------
   if ( u==1 )
   {
      setring R;
      module D2 = fetch(@R,D2);
      if( opt==1 or opt==2 )
      { module D1 = fetch(@R,D1);}
      kill @R;
   }
   if( notInvertible == 1 )
   {
     // The matrix A seems to be non-invertible.
     // Note that there are examples, where this is not true but only due to
     // inexact computations in the field of reals or complex numbers:
     // ring r = complex, x, dp;
     // The following matrix has non-zero determinante but seems non-invertible:
     // matrix A[3][3] = 1,i,i,0,1,2,1,0,1+i;
     // For this example, inverse_B yields the correct answer.
     // So, let's use this as a workaround whenever we have this situation:
     list myList = inverse_B(A);
     matrix Try = inverse_B(A)[1];
     if (myList[2] == poly(1)) { return (Try); }
     else
     {
       dbprint(db,"// ** matrix is not invertible");
       return(invA);
     }
   }
   else { return(matrix(D2)); }     //matrix invertible with inverse D2

}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0,(x,y,z),lp;
  matrix A[3][3]=
   1,4,3,
   1,5,7,
   0,4,17;
  print(inverse(A));"";
  matrix B[3][3]=
   y+1,  x+y,    y,
   z,    z+1,    z,
   y+z+2,x+y+z+2,y+z+1;
  print(inverse(B));
  print(B*inverse(B));
}

//////////////////////////////////////////////////////////////////////////////
proc sym_gauss(matrix A)
"USAGE:    sym_gauss(A);  A = symmetric matrix
RETURN:   matrix, diagonalisation of A with symmetric gauss algorithm
EXAMPLE:  example sym_gauss; shows an example"
{
  int i,j;
  int n=nrows(A);

  if (ncols(A)!=n){
    "// ** input is not a square matrix";;
    return(A);
  }

  if(!const_mat(A)){
    "// ** input is not a constant matrix";
    return(A);
  }

  if(deg(std(A-transpose(A))[1])!=-1){
    "// ** input is not a symmetric matrix";
    return(A);
  }

  for(i=1; i<n; i++){
    for(j=i+1; j<=n; j++){
      if(A[i,j]!=0){ A=red(A,i,j); }
    }
  }

  return(A);
}
example
{"EXAMPLE:"; echo = 2;
  ring r=0,(x),lp;
  matrix A[2][2]=1,4,4,15;
  print(A);
  print(sym_gauss(A));
}

//////////////////////////////////////////////////////////////////////////////
proc orthogonalize(matrix A)
"USAGE:    orthogonalize(A); A = matrix of constants
RETURN:    matrix, orthogonal basis of the colum space of A
EXAMPLE:   example orthogonalize; shows an example "
{
  int i,j;
  int n=ncols(A);
  poly k;

  if(!const_mat(A)){
    "// ** input is not a constant matrix";
    matrix B;
    return(B);
  }

  module B=module(interred(A));

  for(i=1;i<=n;i=i+1) {
    for(j=1;j<i;j=j+1) {
      k=inner_product(B[j],B[j]);
      if (k==0) { "Error: vector of length zero"; return(matrix(B)); }
      B[i]=B[i]-(inner_product(B[i],B[j])/k)*B[j];
    }
  }

  return(matrix(B));
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0,(x),lp;
  matrix A[4][4]=5,6,12,4,7,3,2,6,12,1,1,2,6,4,2,10;
  print(A);
  print(orthogonalize(A));
}

////////////////////////////////////////////////////////////////////////////
proc diag_test(matrix A)
"USAGE:          diag_test(A); A = const square matrix
RETURN:   int,  1 if A is diagonalizable,@*
                0 if not@*
               -1 if no statement is possible, since A does not split.
NOTE:     The test works only for split matrices, i.e if eigenvalues of A
          are in the ground field.
          Does not work with parameters (uses factorize,gcd).
EXAMPLE:  example diag_test; shows an example"
{
  int i,j;
  int n     = nrows(A);
  string mp = string(minpoly);
  string cs = charstr(basering);
  int np=0;

  if(ncols(A) != n) {
    "// input is not a square matrix";
    return(-1);
  }

   if(!const_mat(A)){
    "// input is not a constant matrix";
    return(-1);
  }

  //Parameterring wegen factorize nicht erlaubt
  for(i=1;i<size(cs);i=i+1){
    if(cs[i]==","){np=np+1;} //Anzahl der Parameter
  }
  if(np>0){
        "// rings with parameters not allowed";
        return(-1);
  }

  //speichern des aktuellen Rings
  def BR=basering;
  //setze R[t]
  execute("ring rt=("+charstr(basering)+"),(@t,"+varstr(basering)+"),lp;");
  execute("minpoly="+mp+";");
  matrix A=imap(BR,A);

  intvec z;
  intvec s;
  poly X;         //characteristisches Polynom
  poly dXdt;      //Ableitung von X nach t
  ideal g;              //ggT(X,dXdt)
  poly b;              //Komponente der Busadjunkten-Matrix
  matrix E[n][n]; //Einheitsmatrix

  E=E+1;
  A=E*@t-A;
  X=det(A);

  matrix Xfactors=matrix(factorize(X,1));        //zerfaellt die Matrtix ?
  int nf=ncols(Xfactors);

  for(i=1;i<=nf;i++){
    if(lead(Xfactors[1,i])>=@t^2){
      //" matrix does not split";
      setring BR;
      return(-1);
    }
  }

  dXdt=diff(X,@t);
  g=std(ideal(gcd(X,dXdt)));

  //Busadjunkte
  z=2..n;
  for(i=1;i<=n;i++){
    s=2..n;
    for(j=1;j<=n;j++){
      b=det(submat(A,z,s));

      if(0!=reduce(b,g)){
        //" matrix not diagonalizable";
        setring BR;
        return(0);
      }

      s[j]=j;
    }
    z[i]=i;
  }

  //"Die Matrix ist diagonalisierbar";
  setring BR;
  return(1);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0,(x),dp;
  matrix A[4][4]=6,0,0,0,0,0,6,0,0,6,0,0,0,0,0,6;
  print(A);
  diag_test(A);
}

//////////////////////////////////////////////////////////////////////////////
proc busadj(matrix A)
"USAGE:   busadj(A);  A = square matrix (of size nxn)
RETURN:  list L:
@format
         L[1] contains the (n+1) coefficients of the characteristic
              polynomial X of A, i.e.
              X = L[1][1]+..+L[1][k]*t^(k-1)+..+(L[1][n+1])*t^n
         L[2] contains the n (nxn)-matrices Hk which are the coefficients of
              the busadjoint bA = adjoint(E*t-A) of A, i.e.
              bA = (Hn-1)*t^(n-1)+...+Hk*t^k+...+H0,  ( Hk=L[2][k+1] )
@end format
EXAMPLE: example busadj; shows an example"
{
  int k;
  int n    = nrows(A);
  matrix E = unitmat(n);
  matrix H[n][n];
  matrix B[n][n];
  list bA, X, L;
  poly a;

  if(ncols(A) != n) {
    "input is not a square matrix";
    return(L);
  }

  bA   = E;
  X[1] = 1;
  for(k=1; k<n; k++){
    B  = A*bA[1];              //bA[1] is the last H
    a  = -trace(B)/k;
    H  = B+a*E;
    bA = insert(bA,H);
    X  = insert(X,a);
  }
  B = A*bA[1];
  a = -trace(B)/n;
  X = insert(X,a);

  L = insert(L,bA);
  L = insert(L,X);
  return(L);
}
example
{ "EXAMPLE"; echo = 2;
  ring r = 0,(t,x),lp;
  matrix A[2][2] = 1,x2,x,x2+3x;
  print(A);
  list L = busadj(A);
  poly X = L[1][1]+L[1][2]*t+L[1][3]*t2; X;
  matrix bA[2][2] = L[2][1]+L[2][2]*t;
  print(bA);               //the busadjoint of A;
  print(bA*(t*unitmat(2)-A));
}

//////////////////////////////////////////////////////////////////////////////
proc charpoly(matrix A, list #)
"USAGE:   charpoly(A[,v]); A square matrix, v string, name of a variable
RETURN:  poly, the characteristic polynomial det(E*v-A)
         (default: v=name of last variable)
NOTE:    A must be independent of the variable v. The computation uses det.
         If printlevel>0, det(E*v-A) is diplayed recursively.
EXAMPLE: example charpoly; shows an example"
{
  int n = nrows(A);
  int z = nvars(basering);
  int i,j;
  string v;
  poly X;
  if(ncols(A) != n)
  {
    "// input is not a square matrix";
    return(X);
  }
  //---------------------- test for correct variable -------------------------
  if( size(#)==0 ){
    #[1] = varstr(z);
  }
  if( typeof(#[1]) == "string") { v = #[1]; }
  else
  {
    "// 2nd argument must be a name of a variable not contained in the matrix";
    return(X);
  }
  j=-1;
  for(i=1; i<=z; i++)
  {
    if(varstr(i)==v){j=i;}
  }
  if(j==-1)
  {
    "// "+v+" is not a variable in the basering";
    return(X);
  }
  if ( size(select1(module(A),j)) != 0 )
  {
    "// matrix must not contain the variable "+v;
    "// change to a ring with an extra variable, if necessary."
    return(X);
  }
  //-------------------------- compute charpoly ------------------------------
  X = det(var(j)*unitmat(n)-A);
  if( printlevel-voice+2 >0) { showrecursive(X,var(j));}
  return(X);
}
example
{ "EXAMPLE"; echo=2;
  ring r=0,(x,t),dp;
  matrix A[3][3]=1,x2,x,x2,6,4,x,4,1;
  print(A);
  charpoly(A,"t");
}

//////////////////////////////////////////////////////////////////////////////
proc charpoly_B(matrix A, list #)
"USAGE:   charpoly_B(A[,v]); A square matrix, v string, name of a variable
RETURN:  poly, the characteristic polynomial det(E*v-A)
         (default: v=name of last variable)
NOTE:    A must be constant in the variable v. The computation uses busadj(A).
EXAMPLE: example charpoly_B; shows an example"
{
  int i,j;
  string s,v;
  list L;
  int n     = nrows(A);
  poly X    = 0;
  def BR    = basering;
  string mp = string(minpoly);

  if(ncols(A) != n){
    "// input is not a square matrix";
    return(X);
  }

  //test for correct variable
  if( size(#)==0 ){
    #[1] = varstr(nvars(BR));
  }
  if( typeof(#[1]) == "string"){
     v = #[1];
  }
  else{
    "// 2nd argument must be a name of a variable not contained in the matrix";
    return(X);
  }

  j=-1;
  for(i=1; i<=nvars(BR); i++){
    if(varstr(i)==v){j=i;}
  }
  if(j==-1){
    "// "+v+" is not a variable in the basering";
    return(X);
  }

  //var cannot be in A
  s="Wp(";
  for( i=1; i<=nvars(BR); i++ ){
    if(i!=j){ s=s+"0";}
    else{ s=s+"1";}
    if( i<nvars(BR)) {s=s+",";}
  }
  s=s+")";

  def @R=changeord(s);
  setring @R;
  execute("minpoly="+mp+";");
  matrix A = imap(BR,A);
  for(i=1; i<=n; i++){
    if(deg(lead(A)[i])>=1){
      "// matrix must not contain the variable "+v;
      kill @R;
      setring BR;
      return(X);
    }
  }

  //get coefficients and build the char. poly
  kill @R;
  setring BR;
  L = busadj(A);
  for(i=1; i<=n+1; i++){
    execute("X=X+L[1][i]*"+v+"^"+string(i-1)+";");
  }

  return(X);
}
example
{ "EXAMPLE"; echo=2;
  ring r=0,(x,t),dp;
  matrix A[3][3]=1,x2,x,x2,6,4,x,4,1;
  print(A);
  charpoly_B(A,"t");
}

//////////////////////////////////////////////////////////////////////////////
proc adjoint(matrix A)
"USAGE:    adjoint(A);  A = square matrix
RETURN:   adjoint matrix of A, i.e. Adj*A=det(A)*E
NOTE:     computation uses busadj(A)
EXAMPLE:  example adjoint; shows an example"
{
  int n=nrows(A);
  matrix Adj[n][n];
  list L;

  if(ncols(A) != n) {
    "// input is not a square matrix";
    return(Adj);
  }

  L  = busadj(A);
  Adj= (-1)^(n-1)*L[2][1];
  return(Adj);

}
example
{ "EXAMPLE"; echo=2;
  ring r=0,(t,x),lp;
  matrix A[2][2]=1,x2,x,x2+3x;
  print(A);
  matrix Adj[2][2]=adjoint(A);
  print(Adj);                    //Adj*A=det(A)*E
  print(Adj*A);
}

//////////////////////////////////////////////////////////////////////////////
proc inverse_B(matrix A)
"USAGE:    inverse_B(A);  A = square matrix
RETURN:   list Inv with
          - Inv[1] = matrix I and
          - Inv[2] = poly p
          such that I*A = unitmat(n)*p;
NOTE:     p=1 if 1/det(A) is computable and  p=det(A) if not;
          the computation uses busadj.
SEE ALSO: inverse, inverse_L
EXAMPLE:  example inverse_B; shows an example"
{
  int i;
  int n=nrows(A);
  matrix I[n][n];
  poly factor;
  list L;
  list Inv;

  if(ncols(A) != n) {
    "input is not a square matrix";
    return(I);
  }

  L=busadj(A);
  I=module(-L[2][1]);        //+-Adj(A)

  if(reduce(1,std(L[1][1]))==0){
    I=I*lift(L[1][1],1)[1][1];
    factor=1;
  }
  else{ factor=L[1][1];}     //=+-det(A) or 1
  Inv=insert(Inv,factor);
  Inv=insert(Inv,matrix(I));

  return(Inv);
}
example
{ "EXAMPLE"; echo=2;
  ring r=0,(x,y),lp;
  matrix A[3][3]=x,y,1,1,x2,y,x,6,0;
  print(A);
  list Inv=inverse_B(A);
  print(Inv[1]);
  print(Inv[2]);
  print(Inv[1]*A);
}

//////////////////////////////////////////////////////////////////////////////
proc det_B(matrix A)
"USAGE:     det_B(A);  A any matrix
RETURN:    returns the determinant of A
NOTE:      the computation uses the busadj algorithm
EXAMPLE:   example det_B; shows an example"
{
  int n=nrows(A);
  list L;

  if(ncols(A) != n){ return(0);}

  L=busadj(A);
  return((-1)^n*L[1][1]);
}
example
{ "EXAMPLE"; echo=2;
  ring r=0,(x),dp;
  matrix A[10][10]=random(2,10,10)+unitmat(10)*x;
  print(A);
  det_B(A);
}

//////////////////////////////////////////////////////////////////////////////
proc inverse_L(matrix A)
"USAGE:     inverse_L(A);  A = square matrix
RETURN:    list Inv representing a left inverse of A, i.e
           - Inv[1] = matrix I and
           - Inv[2] = poly p
           such that I*A = unitmat(n)*p;
NOTE:      p=1 if 1/det(A) is computable and p=det(A) if not;
           the computation computes first det(A) and then uses lift
SEE ALSO:  inverse, inverse_B
EXAMPLE:   example inverse_L; shows an example"
{
  int n=nrows(A);
  matrix I;
  matrix E[n][n]=unitmat(n);
  poly factor;
  poly d=1;
  list Inv;

  if (ncols(A)!=n){
    "// input is not a square matrix";
    return(I);
  }

  d=det(A);
  if(d==0){
    "// matrix is not invertible";
    return(Inv);
  }

  // test if 1/det(A) exists
  if(reduce(1,std(d))!=0){ E=E*d;}

  I=lift(A,E);
  if(I==unitmat(n)-unitmat(n)){ //catch error in lift
    "// matrix is not invertible";
    return(Inv);
  }

  factor=d;      //=det(A) or 1
  Inv=insert(Inv,factor);
  Inv=insert(Inv,I);

  return(Inv);
}
example
{ "EXAMPLE"; echo=2;
  ring r=0,(x,y),lp;
  matrix A[3][3]=x,y,1,1,x2,y,x,6,0;
  print(A);
  list Inv=inverse_L(A);
  print(Inv[1]);
  print(Inv[2]);
  print(Inv[1]*A);
}

//////////////////////////////////////////////////////////////////////////////
proc gaussred(matrix A)
"USAGE:   gaussred(A);   A any constant matrix
RETURN:  list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
         gives a row reduced matrix S, a permutation matrix P and a
         normalized lower triangular matrix U, with P*A=U*S
NOTE:    This procedure is designed for teaching purposes mainly.
         The straight forward implementation in the interpreted library
         is not very efficient (no standard basis computation).
EXAMPLE: example gaussred; shows an example"
{
  int i,j,l,k,jp,rang;
  poly c,pivo;
  list Z;
  int n = nrows(A);
  int m = ncols(A);
  int mr= n; //max. rang
  matrix P[n][n] = unitmat(n);
  matrix U[n][n] = P;

  if(!const_mat(A)){
    "// input is not a constant matrix";
    return(Z);
  }

  if(n>m){mr=m;} //max. rang

  for(i=1;i<=mr;i=i+1){
    if((i+k)>m){break;}

    //Test: Diagonalelement=0
    if(A[i,i+k]==0){
      jp=i;pivo=0;
      for(j=i+1;j<=n;j=j+1){
        c=absValue(A[j,i+k]);
        if(pivo<c){ pivo=c;jp=j;}
      }
      if(jp != i){       //Zeilentausch
        for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter)
          c=A[i,j];
          A[i,j]=A[jp,j];
          A[jp,j]=c;
        }
        for(j=1;j<=n;j=j+1){ //Zeilentausch in P
          c=P[i,j];
          P[i,j]=P[jp,j];
          P[jp,j]=c;
        }
      }
      if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe !
    }                          //i sollte im naechsten Lauf nicht erhoeht sein

    //Eliminationsschritt
    for(j=i+1;j<=n;j=j+1){
      c=A[j,i+k]/A[i,i+k];
      for(l=i+k+1;l<=m;l=l+1){
        A[j,l]=A[j,l]-A[i,l]*c;
      }
      A[j,i+k]=0;  // nur wichtig falls k>0 ist
      A[j,i]=c;    // bildet U
    }
  rang=i;
  }

  for(i=1;i<=mr;i=i+1){
    for(j=i+1;j<=n;j=j+1){
      U[j,i]=A[j,i];
      A[j,i]=0;
    }
  }

  Z=insert(Z,rang);
  Z=insert(Z,A);
  Z=insert(Z,U);
  Z=insert(Z,P);

  return(Z);
}
example
{ "EXAMPLE";echo=2;
  ring r=0,(x),dp;
  matrix A[5][4]=1,3,-1,4,2,5,-1,3,1,3,-1,4,0,4,-3,1,-3,1,-5,-2;
  print(A);
  list Z=gaussred(A);   //construct P,U,S s.t. P*A=U*S
  print(Z[1]);          //P
  print(Z[2]);          //U
  print(Z[3]);          //S
  print(Z[4]);          //rank
  print(Z[1]*A);        //P*A
  print(Z[2]*Z[3]);     //U*S
}

//////////////////////////////////////////////////////////////////////////////
proc gaussred_pivot(matrix A)
"USAGE:     gaussred_pivot(A);   A any constant matrix
RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
           gives a row reduced matrix S, a permutation matrix P and a
           normalized lower triangular matrix U, with P*A=U*S
NOTE:      with row pivoting
EXAMPLE:   example gaussred_pivot; shows an example"
{
  int i,j,l,k,jp,rang;
  poly c,pivo;
  list Z;
  int n=nrows(A);
  int m=ncols(A);
  int mr=n; //max. rang
  matrix P[n][n]=unitmat(n);
  matrix U[n][n]=P;

  if(!const_mat(A)){
    "// input is not a constant matrix";
    return(Z);
  }

  if(n>m){mr=m;} //max. rang

  for(i=1;i<=mr;i=i+1){
    if((i+k)>m){break;}

    //Pivotisierung
    pivo=absValue(A[i,i+k]);jp=i;
    for(j=i+1;j<=n;j=j+1){
      c=absValue(A[j,i+k]);
      if(pivo<c){ pivo=c;jp=j;}
    }
    if(jp != i){ //Zeilentausch
      for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter)
        c=A[i,j];
        A[i,j]=A[jp,j];
        A[jp,j]=c;
      }
      for(j=1;j<=n;j=j+1){ //Zeilentausch in P
        c=P[i,j];
        P[i,j]=P[jp,j];
        P[jp,j]=c;
      }
    }
    if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe !
                               //i sollte im naechsten Lauf nicht erhoeht sein
    //Eliminationsschritt
    for(j=i+1;j<=n;j=j+1){
      c=A[j,i+k]/A[i,i+k];
      for(l=i+k+1;l<=m;l=l+1){
        A[j,l]=A[j,l]-A[i,l]*c;
      }
      A[j,i+k]=0;  // nur wichtig falls k>0 ist
      A[j,i]=c;    // bildet U
    }
  rang=i;
  }

  for(i=1;i<=mr;i=i+1){
    for(j=i+1;j<=n;j=j+1){
      U[j,i]=A[j,i];
      A[j,i]=0;
    }
  }

  Z=insert(Z,rang);
  Z=insert(Z,A);
  Z=insert(Z,U);
  Z=insert(Z,P);

  return(Z);
}
example
{ "EXAMPLE";echo=2;
  ring r=0,(x),dp;
  matrix A[5][4] = 1, 3,-1,4,
                   2, 5,-1,3,
                   1, 3,-1,4,
                   0, 4,-3,1,
                  -3,1,-5,-2;
  list Z=gaussred_pivot(A);  //construct P,U,S s.t. P*A=U*S
  print(Z[1]);               //P
  print(Z[2]);               //U
  print(Z[3]);               //S
  print(Z[4]);               //rank
  print(Z[1]*A);             //P*A
  print(Z[2]*Z[3]);          //U*S
}

//////////////////////////////////////////////////////////////////////////////
proc gauss_nf(matrix A)
"USAGE:    gauss_nf(A); A any constant matrix
RETURN:   matrix; gauss normal form of A (uses gaussred)
EXAMPLE:  example gauss_nf; shows an example"
{
  list Z;
  if(!const_mat(A)){
    "// input is not a constant matrix";
    return(A);
  }
  Z = gaussred(A);
  return(Z[3]);
}
example
{ "EXAMPLE";echo=2;
  ring r = 0,(x),dp;
  matrix A[4][4] = 1,4,4,7,2,5,5,4,4,1,1,3,0,2,2,7;
  print(gauss_nf(A));
}

//////////////////////////////////////////////////////////////////////////////
proc mat_rk(matrix A)
"USAGE:    mat_rk(A); A any constant matrix
RETURN:   int, rank of A
EXAMPLE:  example mat_rk; shows an example"
{
  list Z;
  if(!const_mat(A)){
    "// input is not a constant matrix";
    return(-1);
  }
  Z = gaussred(A);
  return(Z[4]);
}
example
{ "EXAMPLE";echo=2;
  ring r = 0,(x),dp;
  matrix A[4][4] = 1,4,4,7,2,5,5,4,4,1,1,3,0,2,2,7;
  mat_rk(A);
}

//////////////////////////////////////////////////////////////////////////////
proc U_D_O(matrix A)
"USAGE:     U_D_O(A);   constant invertible matrix A
RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=D , Z[4]=O
           gives a permutation matrix P,
           a normalized lower triangular matrix U ,
           a diagonal matrix D, and
           a normalized upper triangular matrix O
           with P*A=U*D*O
NOTE:      Z[1]=-1 means that A is not regular (proc uses gaussred)
EXAMPLE:   example U_D_O; shows an example"
{
  int i,j;
  list Z,L;
  int n=nrows(A);
  matrix O[n][n]=unitmat(n);
  matrix D[n][n];

  if (ncols(A)!=n){
    "// input is not a square matrix";
    return(Z);
  }
  if(!const_mat(A)){
    "// input is not a constant matrix";
    return(Z);
  }

  L=gaussred(A);

  if(L[4]!=n){
    "// input is not an invertible matrix";
    Z=insert(Z,-1);  //hint for calling procedures
    return(Z);
  }

  D=L[3];

  for(i=1; i<=n; i++){
    for(j=i+1; j<=n; j++){
      O[i,j] = D[i,j]/D[i,i];
      D[i,j] = 0;
    }
  }

  Z=insert(Z,O);
  Z=insert(Z,D);
  Z=insert(Z,L[2]);
  Z=insert(Z,L[1]);
  return(Z);
}
example
{ "EXAMPLE";echo=2;
  ring r = 0,(x),dp;
  matrix A[5][5] = 10, 4,  0, -9,  8,
                   -3, 6, -6, -4,  9,
                    0, 3, -1, -9, -8,
                   -4,-2, -6, -10,10,
                   -9, 5, -1, -6,  5;
  list Z = U_D_O(A);              //construct P,U,D,O s.t. P*A=U*D*O
  print(Z[1]);                    //P
  print(Z[2]);                    //U
  print(Z[3]);                    //D
  print(Z[4]);                    //O
  print(Z[1]*A);                  //P*A
  print(Z[2]*Z[3]*Z[4]);          //U*D*O
}

//////////////////////////////////////////////////////////////////////////////
proc pos_def(matrix A)
"USAGE:     pos_def(A); A = constant, symmetric square matrix
RETURN:    int:
           1  if A is positive definit ,
           0  if not,
           -1 if unknown
EXAMPLE:   example pos_def; shows an example"
{
  int j;
  list Z;
  int n = nrows(A);
  matrix H[n][n];

  if (ncols(A)!=n){
    "// input is not a square matrix";
    return(0);
  }
  if(!const_mat(A)){
    "// input is not a constant matrix";
    return(-1);
  }
  if(deg(std(A-transpose(A))[1])!=-1){
    "// input is not a hermitian (symmetric) matrix";
    return(-1);
  }

  Z=U_D_O(A);

  if(Z[1]==-1){
    return(0);
  }  //A not regular, therefore not pos. definit

  H=Z[1];
  //es fand Zeilentausch statt: also nicht positiv definit
  if(deg(std(H-unitmat(n))[1])!=-1){
    return(0);
  }

  H=Z[3];

  for(j=1;j<=n;j=j+1){
    if(H[j,j]<=0){
      return(0);
    } //eigenvalue<=0, not pos.definit
  }

  return(1); //positiv definit;
}
example
{ "EXAMPLE"; echo=2;
  ring r = 0,(x),dp;
  matrix A[5][5] = 20,  4,  0, -9,   8,
                    4, 12, -6, -4,   9,
                    0, -6, -2, -9,  -8,
                   -9, -4, -9, -20, 10,
                    8,  9, -8,  10, 10;
  pos_def(A);
  matrix B[3][3] =  3,  2,  0,
                    2, 12,  4,
                    0,  4,  2;
  pos_def(B);
}

//////////////////////////////////////////////////////////////////////////////
proc linsolve(matrix A, matrix b)
"USAGE:     linsolve(A,b); A a constant nxm-matrix, b a constant nx1-matrix
RETURN:    a 1xm matrix X, solution of inhomogeneous linear system A*X = b
           return the 0-matrix if system is not solvable
NOTE:      uses gaussred
EXAMPLE:   example linsolve; shows an example"
{
  int i,j,k,rc,r;
  poly c;
  list Z;
  int n  = nrows(A);
  int m  = ncols(A);
  int n_b= nrows(b);
  matrix Ab[n][m+1];
  matrix X[m][1];

  if(ncols(b)!=1){
    "// right hand side b is not a nx1 matrix";
    return(X);
  }

  if(!const_mat(A)){
    "// input hand is not a constant matrix";
    return(X);
  }

  if(n_b>n){
    for(i=n; i<=n_b; i++){
      if(b[i,1]!=0){
        "// right hand side b not in Image(A)";
        return X;
      }
    }
  }

  if(n_b<n){
    matrix copy[n_b][1]=b;
    matrix b[n][1]=0;
    for(i=1;i<=n_b;i=i+1){
      b[i,1]=copy[i,1];
    }
  }

  r=mat_rk(A);

  //1. b constant vector
  if(const_mat(b)){
    //extend A with b
    for(i=1; i<=n; i++){
      for(j=1; j<=m; j++){
        Ab[i,j]=A[i,j];
      }
      Ab[i,m+1]=b[i,1];
    }

    //Gauss reduction
    Z  = gaussred(Ab);
    Ab = Z[3];  //normal form
    rc = Z[4];  //rank(Ab)
    //print(Ab);

    if(r<rc){
        "// no solution";
        return(X);
    }
    k=m;
    for(i=r;i>=1;i=i-1){

      j=1;
      while(Ab[i,j]==0){j=j+1;}// suche Ecke

      for(;k>j;k=k-1){ X[k]=0;}//springe zur Ecke


      c=Ab[i,m+1]; //i-te Komponene von b
      for(j=m;j>k;j=j-1){
        c=c-X[j,1]*Ab[i,j];
      }
      if(Ab[i,k]==0){
        X[k,1]=1; //willkuerlich
      }
      else{
          X[k,1]=c/Ab[i,k];
      }
      k=k-1;
      if(k==0){break;}
    }


  }//endif (const b)
  else{  //b not constant
    "// !not implemented!";

  }

  return(X);
}
example
{ "EXAMPLE";echo=2;
  ring r=0,(x),dp;
  matrix A[3][2] = -4,-6,
                    2, 3,
                   -5, 7;
  matrix b[3][1] = 10,
                   -5,
                    2;
  matrix X = linsolve(A,b);
  print(X);
  print(A*X);
}
//////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////
//    PROCEDURES for Jordan normal form
//    AUTHOR:  Mathias Schulze, email: mschulze@mathematik.uni-kl.de
///////////////////////////////////////////////////////////////////////////////

static proc rowcolswap(matrix M,int i,int j)
{
  if(i==j)
  {
    return(M);
  }
  poly p;
  for(int k=1;k<=nrows(M);k++)
  {
    p=M[i,k];
    M[i,k]=M[j,k];
    M[j,k]=p;
  }
  for(k=1;k<=ncols(M);k++)
  {
    p=M[k,i];
    M[k,i]=M[k,j];
    M[k,j]=p;
  }
  return(M);
}
//////////////////////////////////////////////////////////////////////////////

static proc rowelim(matrix M,int i,int j,int k)
{
  if(jet(M[i,k],0)==0||jet(M[j,k],0)==0)
  {
    return(M);
  }
  number n=number(jet(M[i,k],0))/number(jet(M[j,k],0));
  for(int l=1;l<=ncols(M);l++)
  {
    M[i,l]=M[i,l]-n*M[j,l];
  }
  for(l=1;l<=nrows(M);l++)
  {
    M[l,j]=M[l,j]+n*M[l,i];
  }
  return(M);
}
///////////////////////////////////////////////////////////////////////////////

static proc colelim(matrix M,int i,int j,int k)
{
  if(jet(M[k,i],0)==0||jet(M[k,j],0)==0)
  {
    return(M);
  }
  number n=number(jet(M[k,i],0))/number(jet(M[k,j],0));
  for(int l=1;l<=nrows(M);l++)
  {
    M[l,i]=M[l,i]-n*M[l,j];
  }
  for(l=1;l<=ncols(M);l++)
  {
    M[j,l]=M[j,l]+n*M[i,l];
  }
  return(M);
}
///////////////////////////////////////////////////////////////////////////////

proc hessenberg(matrix M)
"USAGE:   hessenberg(M); matrix M
ASSUME:  M constant square matrix
RETURN:  matrix H;  Hessenberg form of M
EXAMPLE: example hessenberg; shows examples
"
{
  if(system("with","eigenval"))
  {
    return(system("hessenberg",M));
  }

  int n=ncols(M);
  int i,j;
  for(int k=1;k<n-1;k++)
  {
    j=k+1;
    while(j<n&&jet(M[j,k],0)==0)
    {
      j++;
    }
    if(jet(M[j,k],0)!=0)
    {
      M=rowcolswap(M,j,k+1);
      for(i=j+1;i<=n;i++)
      {
        M=rowelim(M,i,k+1,k);
      }
    }
  }
  return(M);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,x,dp;
  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
  print(M);
  print(hessenberg(M));
}
///////////////////////////////////////////////////////////////////////////////

proc eigenvals(matrix M)
"USAGE:   eigenvals(M); matrix M
ASSUME:  eigenvalues of M in basefield
RETURN:
@format
list l;
  ideal l[1];
    number l[1][i];  i-th eigenvalue of M
  intvec l[2];
    int l[2][i];  multiplicity of i-th eigenvalue of M
@end format
EXAMPLE: example eigenvals; shows examples
"
{
  if(system("with","eigenval"))
  {
    return(system("eigenvals",jet(M,0)));
  }

  M=jet(hessenberg(M),0);
  int n=ncols(M);
  int k;
  ideal e;
  intvec m;
  number e0;
  intvec v;
  list l;
  int i,j;
  j=1;
  while(j<=n)
  {
    v=j;
    j++;
    if(j<=n)
    {
      while(j<n&&M[j,j-1]!=0)
      {
        v=v,j;
        j++;
      }
      if(M[j,j-1]!=0)
      {
        v=v,j;
        j++;
      }
    }
    if(size(v)==1)
    {
      k++;
      e[k]=M[v,v];
      m[k]=1;
    }
    else
    {
      l=factorize(det(submat(M,v,v)-var(1)));
      for(i=size(l[1]);i>=1;i--)
      {
        e0=number(jet(l[1][i]/var(1),0));
        if(e0!=0)
        {
          k++;
          e[k]=(e0*var(1)-l[1][i])/e0;
          m[k]=l[2][i];
        }
      }
    }
  }
  return(spnf(list(e,m)));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,x,dp;
  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
  print(M);
  eigenvals(M);
}
///////////////////////////////////////////////////////////////////////////////

proc minipoly(matrix M,list #)
"USAGE:   minipoly(M); matrix M
ASSUME:  eigenvalues of M in basefield
RETURN:
@format
list l;  minimal polynomial of M
  ideal l[1];
    number l[1][i];  i-th root of minimal polynomial of M
  intvec l[2];
    int l[2][i];  multiplicity of i-th root of minimal polynomial of M
@end format
EXAMPLE: example minipoly; shows examples
"
{
  if(nrows(M)==0)
  {
    ERROR("non empty expected");
  }
  if(ncols(M)!=nrows(M))
  {
    ERROR("square matrix expected");
  }

  M=jet(M,0);

  if(size(#)==0)
  {
    #=eigenvals(M);
  }
  def e0,m0=#[1..2];

  intvec m1;
  matrix N0,N1;
  for(int i=1;i<=ncols(e0);i++)
  {
    m1[i]=1;
    N0=M-e0[i];
    N1=N0;
    while(size(syz(N1))<m0[i])
    {
      m1[i]=m1[i]+1;
      N1=N1*N0;
    }
  }

  return(list(e0,m1));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,x,dp;
  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
  print(M);
  minipoly(M);
}
///////////////////////////////////////////////////////////////////////////////

proc spnf(list #)
"USAGE:   spnf(list(a[,m])); ideal a, intvec m
ASSUME:  ncols(a)==size(m)
RETURN:  list l:
            l[1] an ideal, the generators of a; sorted and with multiple entries displayed only once@*
            l[2] and intvec, l[2][i] provides the multiplicity of l[1][i]
EXAMPLE: example spnf; shows examples
"
{
  list sp=#;
  ideal a=sp[1];
  int n=ncols(a);
  intvec m;
  list V;
  module v;
  int i,j;
  for(i=2;i<=size(sp);i++)
  {
    if(typeof(sp[i])=="intvec")
    {
      m=sp[i];
    }
    if(typeof(sp[i])=="module")
    {
      v=sp[i];
      for(j=n;j>=1;j--)
      {
        V[j]=module(v[j]);
      }
    }
    if(typeof(sp[i])=="list")
    {
      V=sp[i];
    }
  }
  if(m==0)
  {
    for(i=n;i>=1;i--)
    {
      m[i]=1;
    }
  }

  int k;
  ideal a0;
  intvec m0;
  list V0;
  number a1;
  int m1;
  for(i=n;i>=1;i--)
  {
    if(m[i]!=0)
    {
      for(j=i-1;j>=1;j--)
      {
        if(m[j]!=0)
        {
          if(number(a[i])>number(a[j]))
          {
            a1=number(a[i]);
            a[i]=a[j];
            a[j]=a1;
            m1=m[i];
            m[i]=m[j];
            m[j]=m1;
            if(size(V)>0)
            {
              v=V[i];
              V[i]=V[j];
              V[j]=v;
            }
          }
          if(number(a[i])==number(a[j]))
          {
            m[i]=m[i]+m[j];
            m[j]=0;
            if(size(V)>0)
            {
              V[i]=V[i]+V[j];
            }
          }
        }
      }
      k++;
      a0[k]=a[i];
      m0[k]=m[i];
      if(size(V)>0)
      {
        V0[k]=V[i];
      }
    }
  }

  if(size(V0)>0)
  {
    n=size(V0);
    module U=std(V0[n]);
    for(i=n-1;i>=1;i--)
    {
      V0[i]=simplify(reduce(V0[i],U),1);
      if(i>=2)
      {
        U=std(U+V0[i]);
      }
    }
  }

  if(k>0)
  {
    sp=a0,m0;
    if(size(V0)>0)
    {
      sp[3]=V0;
    }
  }
  return(sp);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list sp=list(ideal(-1/2,-3/10,-3/10,-1/10,-1/10,0,1/10,1/10,3/10,3/10,1/2));
  spprint(spnf(sp));
}
///////////////////////////////////////////////////////////////////////////////

proc spprint(list sp)
"USAGE:   spprint(sp); list sp (helper routine for spnf)
RETURN:  string s;  spectrum sp
EXAMPLE: example spprint; shows examples
SEE ALSO: gmssing_lib, spnf
"
{
  string s;
  for(int i=1;i<size(sp[2]);i++)
  {
    s=s+"("+string(sp[1][i])+","+string(sp[2][i])+"),";
  }
  s=s+"("+string(sp[1][i])+","+string(sp[2][i])+")";
  return(s);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,(x,y),ds;
  list sp=list(ideal(-1/2,-3/10,-1/10,0,1/10,3/10,1/2),intvec(1,2,2,1,2,2,1));
  spprint(sp);
}
///////////////////////////////////////////////////////////////////////////////

proc jordan(matrix M,list #)
"USAGE:   jordan(M); matrix M
ASSUME:  eigenvalues of M in basefield
RETURN:
@format
list l;  Jordan data of M
  ideal l[1];
    number l[1][i];  eigenvalue of i-th Jordan block of M
  intvec l[2];
    int l[2][i];  size of i-th Jordan block of M
  intvec l[3];
    int l[3][i];  multiplicity of i-th Jordan block of M
@end format
EXAMPLE: example jordan; shows examples
"
{
  if(nrows(M)==0)
  {
    ERROR("non empty expected");
  }
  if(ncols(M)!=nrows(M))
  {
    ERROR("square matrix expected");
  }

  M=jet(M,0);

  if(size(#)==0)
  {
    #=eigenvals(M);
  }
  def e0,m0=#[1..2];

  int i;
  for(i=1;i<=ncols(e0);i++)
  {
    if(deg(e0[i])>0)
    {

      ERROR("eigenvalues in coefficient field expected");
      return(list());
    }
  }

  int j,k;
  matrix N0,N1;
  module K0;
  list K;
  ideal e;
  intvec s,m;

  for(i=1;i<=ncols(e0);i++)
  {
    N0=M-e0[i]*matrix(freemodule(ncols(M)));

    N1=N0;
    K0=0;
    K=module();
    while(size(K0)<m0[i])
    {
      K0=syz(N1);
      K=K+list(K0);
      N1=N1*N0;
    }

    for(j=2;j<size(K);j++)
    {
      if(2*size(K[j])-size(K[j-1])-size(K[j+1])>0)
      {
        k++;
        e[k]=e0[i];
        s[k]=j-1;
        m[k]=2*size(K[j])-size(K[j-1])-size(K[j+1]);
      }
    }
    if(size(K[j])-size(K[j-1])>0)
    {
      k++;
      e[k]=e0[i];
      s[k]=j-1;
      m[k]=size(K[j])-size(K[j-1]);
    }
  }

  return(list(e,s,m));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,x,dp;
  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
  print(M);
  jordan(M);
}
///////////////////////////////////////////////////////////////////////////////

proc jordanbasis(matrix M,list #)
"USAGE:   jordanbasis(M); matrix M
ASSUME:  eigenvalues of M in basefield
RETURN:
@format
list l:
  module l[1];  inverse(l[1])*M*l[1] in Jordan normal form
  intvec l[2];
    int l[2][i];  weight filtration index of l[1][i]
@end format
EXAMPLE: example jordanbasis; shows examples
"
{
  if(nrows(M)==0)
  {
    ERROR("non empty matrix expected");
  }
  if(ncols(M)!=nrows(M))
  {
    ERROR("square matrix expected");
  }

  M=jet(M,0);

  if(size(#)==0)
  {
    #=eigenvals(M);
  }
  def e,m=#[1..2];

  for(int i=1;i<=ncols(e);i++)
  {
    if(deg(e[i])>0)
    {
      ERROR("eigenvalues in coefficient field expected");
      return(freemodule(ncols(M)));
    }
  }

  int j,k,l,n;
  matrix N0,N1;
  module K0,K1;
  list K;
  matrix u[ncols(M)][1];
  module U;
  intvec w;

  for(i=1;i<=ncols(e);i++)
  {
    N0=M-e[i]*matrix(freemodule(ncols(M)));

    N1=N0;
    K0=0;
    K=list();
    while(size(K0)<m[i])
    {
      K0=syz(N1);
      K=K+list(K0);
      N1=N1*N0;
    }

    K1=0;
    for(j=1;j<size(K);j++)
    {
      K0=K[j];
      K[j]=interred(reduce(K[j],std(K1+module(N0*K[j+1]))));
      K1=K0;
    }
    K[j]=interred(reduce(K[j],std(K1)));

    for(l=size(K);l>=1;l--)
    {
      for(k=size(K[l]);k>0;k--)
      {
        u=K[l][k];
        for(j=l;j>=1;j--)
        {
          U=U+module(u);
          n++;
          w[n]=2*j-l-1;
          u=N0*u;
        }
      }
    }
  }

  return(list(U,w));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,x,dp;
  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
  print(M);
  list l=jordanbasis(M);
  print(l[1]);
  print(l[2]);
  print(inverse(l[1])*M*l[1]);
}
///////////////////////////////////////////////////////////////////////////////

proc jordanmatrix(list jd)
"USAGE:   jordanmatrix(list(e,s,m)); ideal e, intvec s, intvec m
ASSUME:  ncols(e)==size(s)==size(m)
RETURN:
@format
matrix J;  Jordan matrix with list(e,s,m)==jordan(J)
@end format
EXAMPLE: example jordanmatrix; shows examples
"
{
  ideal e=jd[1];
  intvec s=jd[2];
  intvec m=jd[3];
  if(ncols(e)!=size(s)||ncols(e)!=size(m))
  {
    ERROR("arguments of equal size expected");
  }

  int i,j,k,l;
  int n=int((transpose(matrix(s))*matrix(m))[1,1]);
  matrix J[n][n];
  for(k=1;k<=ncols(e);k++)
  {
    for(l=1;l<=m[k];l++)
    {
      j++;
      J[j,j]=e[k];
      for(i=s[k];i>=2;i--)
      {
        J[j+1,j]=1;
        j++;
        J[j,j]=e[k];
      }
    }
  }

  return(J);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,x,dp;
  ideal e=ideal(2,3);
  intvec s=1,2;
  intvec m=1,1;
  print(jordanmatrix(list(e,s,m)));
}
///////////////////////////////////////////////////////////////////////////////

proc jordannf(matrix M,list #)
"USAGE:   jordannf(M); matrix M
ASSUME:  eigenvalues of M in basefield
RETURN:  matrix J; Jordan normal form of M
EXAMPLE: example jordannf; shows examples
"
{
  return(jordanmatrix(jordan(M,#)));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,x,dp;
  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
  print(M);
  print(jordannf(M));
}

///////////////////////////////////////////////////////////////////////////////

/*
///////////////////////////////////////////////////////////////////////////////
//          Auskommentierte zusaetzliche Beispiele
//
///////////////////////////////////////////////////////////////////////////////
// Singular for ix86-Linux version 1-3-10  (2000121517)  Dec 15 2000 17:55:12
// Rechnungen auf AMD700 mit 632 MB

 LIB "linalg.lib";

1. Sparse integer Matrizen
--------------------------
ring r1=0,(x),dp;
system("--random", 12345678);
int n = 70;
matrix m = sparsemat(n,n,50,100);
option(prot,mem);

int t=timer;
matrix im = inverse(m,1)[1];
timer-t;
print(im*m);
//list l0 = watchdog(100,"inverse("+"m"+",3)");
//bricht bei 100 sec ab und gibt l0[1]: string Killed zurueck

//inverse(m,1): std       5sec 5,5 MB
//inverse(m,2): interred  12sec
//inverse(m,2): lift      nach 180 sec 13MB abgebrochen
//n=60: linalgorig: 3  linalg: 5
//n=70: linalgorig: 6,7 linalg: 11,12
// aber linalgorig rechnet falsch!

2. Sparse poly Matrizen
-----------------------
ring r=(0),(a,b,c),dp;
system("--random", 12345678);
int n=6;
matrix m = sparsematrix(n,n,2,0,50,50,9); //matrix of polys of deg <=2
option(prot,mem);

int t=timer;
matrix im = inverse(m);
timer-t;
print(im*m);
//inverse(m,1): std       0sec 1MB
//inverse(m,2): interred  0sec 1MB
//inverse(m,2): lift      nach 2000 sec 33MB abgebrochen

3. Sparse Matrizen mit Parametern
---------------------------------
//liborig rechnet hier falsch!
ring r=(0),(a,b),dp;
system("--random", 12345678);
int n=7;
matrix m = sparsematrix(n,n,1,0,40,50,9);
ring r1 = (0,a,b),(x),dp;
matrix m = imap(r,m);
option(prot,mem);

int t=timer;
matrix im = inverse(m);
timer-t;
print(im*m);
//inverse(m)=inverse(m,3):15 sec inverse(m,1)=1sec inverse(m,2):>120sec
//Bei Parametern vergeht die Zeit beim Normieren!

3. Sparse Matrizen mit Variablen und Parametern
-----------------------------------------------
ring r=(0),(a,b),dp;
system("--random", 12345678);
int n=6;
matrix m = sparsematrix(n,n,1,0,35,50,9);
ring r1 = (0,a),(b),dp;
matrix m = imap(r,m);
option(prot,mem);

int t=timer;
matrix im = inverse(m,3);
timer-t;
print(im*m);
//n=7: inverse(m,3):lange sec inverse(m,1)=1sec inverse(m,2):1sec

4. Ueber Polynomring invertierbare Matrizen
-------------------------------------------
LIB"random.lib"; LIB"linalg.lib";
system("--random", 12345678);
int n =3;
ring r= 0,(x,y,z),(C,dp);
matrix A=triagmatrix(n,n,1,0,0,50,2);
intmat B=sparsetriag(n,n,20,1);
matrix M = A*transpose(B);
M=M*transpose(M);
M[1,1..ncols(M)]=M[1,1..n]+xyz*M[n,1..ncols(M)];
print(M);
//M hat det=1 nach Konstruktion

int t=timer;
matrix iM=inverse(M);
timer-t;
print(iM*M);                      //test

//ACHTUNG: Interred liefert i.A. keine Inverse, Gegenbeispiel z.B.
//mit n=3
//eifacheres Gegenbeispiel:
matrix M =
9yz+3y+3z+2,             9y2+6y+1,
9xyz+3xy+3xz-9z2+2x-6z-1,9xy2+6xy-9yz+x-3y-3z
//det M=1, inverse(M,2); ->// ** matrix is not invertible
//lead(M); 9xyz*gen(2) 9xy2*gen(2) nicht teilbar!

5. charpoly:
-----------
//ring rp=(0,A,B,C),(x),dp;
ring r=0,(A,B,C,x),dp;
matrix m[12][12]=
AC,BC,-3BC,0,-A2+B2,-3AC+1,B2,   B2,  1,   0, -C2+1,0,
1, 1, 2C,  0,0,     B,     -A,   -4C, 2A+1,0, 0,    0,
0, 0, 0,   1,0,     2C+1,  -4C+1,-A,  B+1, 0, B+1,  3B,
AB,B2,0,   1,0,     1,     0,    1,   A,   0, 1,    B+1,
1, 0, 1,   0,0,     1,     0,    -C2, 0,   1, 0,    1,
0, 0, 2,   1,2A,    1,     0,    0,   0,   0, 1,    1,
0, 1, 0,   1,1,     2,     A,    3B+1,1,   B2,1,    1,
0, 1, 0,   1,1,     1,     1,    1,   2,   0, 0,    0,
1, 0, 1,   0,0,     0,     1,    0,   1,   1, 0,    3,
1, 3B,B2+1,0,0,     1,     0,    1,   0,   0, 1,    0,
0, 0, 1,   0,0,     0,     0,    1,   0,   0, 0,    0,
0, 1, 0,   1,1,     3,     3B+1, 0,   1,   1, 1,    0;
option(prot,mem);

int t=timer;
poly q=charpoly(m,"x");         //1sec, charpoly_B 1sec, 16MB
timer-t;
//1sec, charpoly_B 1sec, 16MB  (gleich in r und rp)

*/