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

version="$Id: jordan.lib,v 1.11 1999-06-16 17:24:35 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, Jordan transformation 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
 invmat(M);        inverse matrix of invertible constant matrix M
";

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

static proc countblocks(matrix M)
{
  int n;

  int i=1;
  int r,r0;
  while(i<=nrows(M))
  {
    n++;
    r=nrows(M[i]);
    r0=r;

    while(i<r0&&i<nrows(M))
    {
      i++;
      if(i<=nrows(M))
      {
        r=nrows(M[i]);
        if(r>r0)
        {
          r0=r;
        }
      }
    }

    i++;
  }

  return(n);
}
///////////////////////////////////////////////////////////////////////////////

static proc getblock(matrix M,intvec v)
{
  matrix M0[size(v)][size(v)]=M[v,v];
  return(M0);
}
///////////////////////////////////////////////////////////////////////////////

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:  The procedure returns a list jd with 3 entries of type
         ideal, list of intvecs, matrix with
         jd[1] eigenvalues of M,
         jd[2][i][j] size of j-th Jordan block with eigenvalue jd[1][i], and
         jd[3]^(-1)*M*jd[3] in Jordan normal form.
         Depending on opt, only certain entries of jd are computed.
         If opt=-1, jd[1] is computed,
         if opt= 0, jd[1] and jd[2] are computed,
         if opt= 1, jd[1],jd[2], and jd[3] are computed, and,
         if opt= 2, jd[1] and jd[3] are computed.
         By default, opt=0.
NOTE:    A non constant polynomial matrix M is replaced by its constant part.
EXAMPLE: example jordan; shows an example.
"
{
  // check parameter M
  int n=nrows(M);
  if(n==0)
  {
    print("//empty matrix");
    return(list());
  }
  if(n!=ncols(M))
  {
    print("//no square matrix");
    return(list());
  }

  // replace M by its constant part
  M=jet(M,0);

  // try to increase number of blocks by transposing M
  if(countblocks(M)<countblocks(transpose(M)))
  {
    M=transpose(M);
  }

  // compute eigenvalues eM and multiplicities mM of M
  list fac;
  matrix M0;
  ideal eM,eM0;
  intvec mM,mM0;
  intvec u;
  int i=1;
  int j,r,r0;
  while(i<=nrows(M))
  {
    u=i;
    r=nrows(M[i]);
    r0=r;

    // find next block
    while(i<r0&&i<nrows(M))
    {
      i++;
      if(i<=nrows(M))
      {
        u=u,i;
        r=nrows(M[i]);
        if(r>r0)
        {
          r0=r;
        }
      }
    }

    // if 1x1-block
    if(size(u)==1)
    {
      if(mM[1]==0)
      {
        eM=M[u[1]][u[1]];
        mM=1;
      }
      else
      {
        eM=eM,ideal(M[u[1]][u[1]]);
        mM=mM,1;
      }
    }
    else
    {
      // factorize characteristic polynomial of block M0
      M0=getblock(M,u);
      fac=factorize(det(M0-var(1)*freemodule(size(u))),2);
      eM0,mM0=fac[1..2];

      // compute eigenvalues eM0 of M0
      if(1<var(1))
      {
        for(j=ncols(eM0);j>=1;j--)
	{
          if(deg(eM0[j])>1)
          {
            print("//eigenvalues not in the coefficient field");
            return(list());
          }
          eM0[j]=-(eM0[j][2]/var(1))/eM0[j][1];
        }
      }
      else
      {
        for(j=ncols(eM0);j>=1;j--)
	{
          if(deg(eM0[j])>1)
          {
            print("//eigenvalues not in the coefficient field");
            return(list());
          }
          eM0[j]=-eM0[j][1]/(eM0[j][2]/var(1));
        }
      }

      if(mM[1]==0)
      {
        eM=eM0;
        mM=mM0;
      }
      else
      {
        eM=eM,eM0;
        mM=mM,mM0;
      }
    }

    i++;
  }

  // sort eigenvalues
  poly e;
  int m;
  for(i=ncols(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;
      }
    }
  }

  // remove multiple eigenvalues
  i=1;
  j=2;
  while(j<=ncols(eM))
  {
    if(eM[i]==eM[j])
    {
      mM[i]=mM[i]+mM[j];
    }
    else
    {
      i++;
      eM[i]=eM[j];
    }
    j++;
  }
  eM=eM[1..i];
  mM=mM[1..i];

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

  // define 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;
  }

  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 bMi 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 Jordan basis vectors K 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)));
      }
      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;
          }
        }
      }
    }
  }

  list jd=eM;
  if(opt<=1)
  {
    jd[2]=bM;
  }
  if(opt>=1)
  {
    jd[3]=V;
  }
  return(jd);
}
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 jd)
"USAGE:   jordanmatrix(jd); jd list of ideal and list of intvecs
RETURN:  The procedure returns the Jordan matrix J with eigenvalues jd[1] and
         size jd[2][i][j] of j-th Jordan block with eigenvalue jd[1][i].
EXAMPLE: example jordanmatrix; shows an example.
"
{
  // check parameters
  if(size(jd)<2)
  {
    print("//not enough entries in argument list");
    matrix J[1][0];
    return(J); 
  }
  def eJ,bJ=jd[1..2];
  if(typeof(eJ)!="ideal")
  {
    print("//first entry in argument list not an ideal");
    matrix J[1][0];
    return(J); 
  }
  if(typeof(bJ)!="list")
  {
    print("//second entry in argument list not a list");
    matrix J[1][0];
    return(J); 
  }
  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")
    {
      print("//second entry in argument list not a list of intvecs");
      matrix J[1][0];
      return(J);
    }
    else
    {
      for(j=size(bJ[i]);j>=1;j--)
      {
        k=bJ[i][j];
        if(k>0)
        {
          n=n+k;
        }
      }
    }
  }

  // create Jordan matrix
  int l;
  matrix J[n][n];
  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 procedure returns 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 invmat(matrix M)
"USAGE:   invmat(M); M constant square matrix
ASSUME:  M is invertible.
RETURN:  The procedure returns the inverse matrix of M.
NOTE:    A non constant polynomial matrix M is replaced by its constant part.
EXAMPLE: example invmat; shows an example.
"
{
  if(nrows(M)==ncols(M))
  {
    matrix invM=lift(jet(M,0),freemodule(nrows(M)));
  }
  else
  {
    print("//no square matrix");
    matrix[1][0]=invM;
  }
  return(invM);
}
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(invmat(M));
}
///////////////////////////////////////////////////////////////////////////////