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

version="$Id: jordan.lib,v 1.9 1999-02-18 19:29:01 mschulze Exp $";
info="
LIBRARY: jordan.lib  PROCEDURES TO COMPUTE THE JORDAN NORMAL FORM
AUTHOR:  Mathias Schulze, email: mschulze@mathematik.uni-kl.de

 jordan(M[,opt]);  eigenvalues, Jordan block sizes, transformation matrix of M
 jordanmatrix(l);  Jordan matrix with eigenvalues l[1], Jordan block sizes l[2]
 jordanform(M);    Jordan normal form of constant square matrix M
 inversemat(M);    inverse matrix of invertible constant matrix M
";

LIB "ring.lib";
///////////////////////////////////////////////////////////////////////////////

proc jordan(matrix M,list #)
"USAGE:   jordan(M[,opt]); M constant square matrix, opt integer
ASSUME:  the eigenvalues of M are in the coefficient field
RETURN:  a list l with 3 entries of type ideal, list of intvecs, matrix with
           l[1]       : eigenvalues of M
           l[2][i][j] : size of j-th Jordan block with eigenvalue l[1][i]
           l[3] : transformation matrix U with U^(-1)*M*U in Jordan normal form
         which entries of l are computed depends on opt:
           opt=-1 : compute l[1]
           opt= 0 : compute l[1] and l[2] (default)
           opt= 1 : compute l[1],l[2] and l[3]
           opt= 2 : compute l[1] and l[3]
NOTE:    a non constant polynomial matrix M is replaced by its constant part
EXAMPLE: example jordan; shows an example
"
{
  // test if square matrix
  int n=nrows(M);
  if(n!=ncols(M))
  {
    "//no square matrix";
    return();
  }

  // get constant part
  def br=basering;
  map zero=br,0;
  M=zero(M);
  kill zero;

  // change to polynomial ring for factorization
  changeord("pr","dp");
  matrix M=imap(br,M);

  // factorize characteristic polynomial
  list l=factorize(det(M-var(1)*freemodule(n)),2);

  // get multiplicities mM
  def eM,mM=l[1..2];
  kill l;

  // test if eigenvalues in the coefficient field
  int i;
  for(i=size(eM);i>=1;i--)
  {
    if(deg(eM[i])>1)
    {
      kill pr;
      "//eigenvalues not in the coefficient field";
      return();
    }
  }

  // get eigenvalues eM
  map inv=pr,-var(1);
  eM=simplify(inv(eM),1);
  setring br;
  map zero=pr,0;
  ideal eM=zero(eM);
  kill pr;

  // sort eigenvalues
  int j;
  poly e;
  int m;
  for(i=size(eM);i>=2;i--)
  {
    for(j=i-1;j>=1;j--)
    {
      if(eM[i]<eM[j])
      {
        e=eM[i];
        eM[i]=eM[j];
        eM[j]=e;
        m=mM[i];
        mM[i]=mM[j];
        mM[j]=m;
      }
    }
  }
  kill e,m;

  // get option parameter
  int opt=0;
  if(size(#)>0)
  {
    if(typeof(#[1])=="int")
    {
      opt=#[1];
    }
  }
  if(opt<0)
  {
    return(list(eM));
  }

  // define needed variables
  int k,l;
  matrix E=freemodule(n);
  matrix Mi,Ni;
  module sNi;
  list K;
  if(opt>=1)
  {
    module V,K1,K2;
    matrix v[n][1];
  }
  if(opt<=1)
  {
    list bM;
    intvec bMi;
  }

  // do the following for all eigenvalues eM[i]
  for(i=ncols(eM);i>=1;i--)
  {
    Mi=M-eM[i]*E;

    // compute kernels K of powers of Mi
    K=list(module());
    for(Ni,sNi=Mi,0;size(sNi)<mM[i];Ni=Ni*Mi)
    {
      sNi=syz(Ni);
      K=K+list(sNi);
    }

    if(opt<=1)
    {
      // compute Jordan block size vector corresponding to eigenvalue eM[i]
      bMi=0;
      bMi[size(K[2])]=0;
      for(j=size(K);j>=2;j--)
      {
        for(k=size(bMi);k>size(bMi)+size(K[j-1])-size(K[j]);k--)
        {
          bMi[k]=bMi[k]+1;
        }
      }
      bM=list(bMi)+bM;
    }

    if(opt>=1)
    {
      // compute generating vectors for Jordan basis vectors corresponding to 
      // eigenvalue eM[i]
      if(size(K)>1)
      {
        for(j,K1=2,0;j<=size(K)-1;j++)
        {
          K2=K[j];
          K[j]=interred(reduce(K[j],std(K1+module(Mi*K[j+1]))));
          K1=K2;
        }
        K[j]=interred(reduce(K[j],std(K1)));
      }

      // compute Jordan basis vectors corresponding to eigenvalue eM[i] from
      // generating vectors
      for(j=size(K);j>=2;j--)
      {
        for(k=size(K[j]);k>=1;k--)
        {
          v=K[j][k];
          for(l=j;l>=1;l--)
          {
            V=module(v)+V;
            v=Mi*v;
          }
        }
      }
    }
  }
  kill l;

  // create return list
  list l=eM;
  if(opt<=1)
  {
    l[2]=bM;
  }
  if(opt>=1)
  {
    l[3]=V;
  }

  return(l);
}
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 jordanmatrix(list l)
"USAGE:   jordanmatrix(l); l list of ideal and list of intvecs
RETURN:  the Jordan matrix J defined by l:
           l[1]       : eigenvalues of J
           l[2][i][j] : size of j-th Jordan block with eigenvalue l[1][i]
EXAMPLE: example jordanmatrix; shows an example
"
{
  // check parameters
  if(size(l)<2)
  {
    "//not enough entries in argument list";
    return(); 
  }
  def eJ,bJ=l[1..2];
  kill l;
  if(typeof(eJ)!="ideal")
  {
    "//first entry in argument list not an ideal";
    return();
  }
  if(typeof(bJ)!="list")
  {
    "//second entry in argument list not a list";
    return();
  }
  if(size(eJ)<size(bJ))
  {
    int s=size(eJ);
  }
  else
  {
    int s=size(bJ);
  }

  // get size of Jordan matrix
  int i,j,k,n;
  for(i=s;i>=1;i--)
  {
    if(typeof(bJ[i])!="intvec")
    {
      "//second entry in argument list not a list of integer vectors";
      return();
    }
    else
    {
      for(j=size(bJ[i]);j>=1;j--)
      {
        k=bJ[i][j];
        if(k>0)
        {
          n=n+k;
        }
      }
    }
  }
  matrix J[n][n];

  // create Jordan matrix
  int l;
  for(i,l=1,1;i<=s;i++)
  {
    for(j=1;j<=size(bJ[i]);j++)
    {
      k=bJ[i][j];
      if(k>0)
      {
        while(k>=2)
        {
          J[l,l]=eJ[i];
          J[l,l+1]=1;
          k,l=k-1,l+1;
        }
        J[l,l]=eJ[i];
        l++;
      }
    }
  }

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

proc jordanform(matrix M)
"USAGE:   jordanform(M); M constant square matrix
ASSUME:  the eigenvalues of M are in the coefficient field
RETURN:  the Jordan normal form of M
NOTE:    a non constant polynomial matrix M is replaced by its constant part
EXAMPLE: example jordanform; shows an example
"
{
  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(jordanform(M));
}
///////////////////////////////////////////////////////////////////////////////

proc inversemat(matrix M)
"USAGE:   inversemat(M); M constant square matrix
ASSUME:  M is invertible
RETURN:  the inverse matrix of M
NOTE:    a non constant polynomial matrix M is replaced by its constant part
EXAMPLE: example inversemat; shows an example
"
{
  // test if square matrix
  int n=nrows(M);
  if(n!=ncols(M))
  {
    "//no square matrix";
    return();
  }

  // get constant part
  def br=basering;
  map zero=br,0;
  M=zero(M);

  // compute inverse
  return(lift(M,freemodule(n)));
}
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(inversemat(M));
}
///////////////////////////////////////////////////////////////////////////////