Changeset dd08428 in git for Singular/LIB

Timestamp:

Dec 23, 1998, 2:13:09 PM (25 years ago)

Author:

Mathias Schulze <mschulze@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

66968388c61cc598cbf981f8917ef5bfcc506cfa

Parents:

7614f91f1e407b440f3a529647a08b4309a850e6

Message:

*mschulze: added procedure to compute the tranformation matrix to jordan normal form


git-svn-id: file:///usr/local/Singular/svn/trunk@2780 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• Singular/LIB/jordan.lib (modified) (7 diffs)

Singular/LIB/jordan.lib

@@ -1 +1 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-version="$Id: jordan.lib,v 1.1 1998-12-21 16:21:05 mschulze Exp $";
+version="$Id: jordan.lib,v 1.2 1998-12-23 13:13:09 mschulze Exp $";
 info="
-LIBRARY: jordan.lib  A PROCEDURE TO COMPUTE THE JORDAN NORMAL FORM
+LIBRARY: jordan.lib  PROCEDURES TO COMPUTE THE JORDAN NORMAL FORM
                      by Mathias Schulze
                      email: mschulze@mathematik.uni-kl.de
 
- jordan(M);          jordan normal form of the constant part of the matrix M
+ jordandata(M);      data of the jordan normal form of the constant part 
+                     of the matrix M
+ jordanmatrix(d);    jordan matrix defined by the data d
+ jordanform(M);      jordan normal form of the constant part of the matrix M
+ jordanbasis(M);     basis with respect to which the constant part 
+                     of the matrix M is in jordan normal form
 ";
 
@@ -13 +18 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-proc jordan(matrix M,list #)
-INPUT:   matrix M[, int i]
-ASSUME:  M square matrix with factorizable characteristic polynomial
+proc jordandata(matrix M)
+INPUT:   square matrix M with factorizable characteristic polynomial 
+         of the constant part M0
 OUTPUT:  nothing, if the assumptions are not fulfilled, else,
-         if i==1, 
-           a list with first entry an ideal containing the eigenvalues 
-           of the constant part of M and second entry a list of decreasingly 
-           ordered int vectors containing the sizes of the corresponding 
-           jordan blocks, 
-         else, 
-           a jordan matrix of the constant part of M
-EXAMPLE: example jordan; shows an example
+         a list with first entry an ideal containing the eigenvalues of M0 
+         in increasing order and second entry a list of increasingly ordered 
+         int vectors containing the corresponding jordan block sizes
+EXAMPLE: example jordandata; shows an example
 {
   int n=nrows(M);
@@ -61 +62 @@
   kill pr;
 
-  int j,k,l,m;
+  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;
+
+  int k,l,m;
   intvec bMi;
   list bM;
@@ -73 +94 @@
     m=k;
     bMi=0;
+    bMi[m]=0;
     for(j=k;j>=1;j--)
     {
@@ -81 +103 @@
       Mj=Mj*Mi;
       k=ncols(syz(Mj));
-      for(l=k-m;l>=1;l--)
+      for(l=size(bMi);l>size(bMi)+m-k;l--)
       {
         bMi[l]=bMi[l]+1;
@@ -90 +112 @@
   }
 
-  if(size(#)>0)
-  {
-    if(#[1]==1)
-    {
-      return(list(eM,bM));
-    }
-  }
-
-  for(i,l,M=size(eM),1,0;i>=1;i--)
+  return(list(eM,bM));
+}
+example
+{
+  "EXAMPLE:";
+  echo=2;
+  LIB "random.lib";
+  ring r;
+  list d=ideal(2,3),list(intvec(1,1,2,2,2),intvec(3,4,5));
+  matrix M=jordanmatrix(d);
+  print(M);
+  int n=nrows(M);
+  matrix U;
+  while(det(U)==0)
+  {
+    U=randommat(n,n,ideal(1)); 
+  }
+  matrix invU=lift(U,freemodule(n));
+  M=invU*M*U;
+  print(M);
+  jordandata(M);
+}
+///////////////////////////////////////////////////////////////////////////////
+
+proc jordanmatrix(list d)
+INPUT:   list d
+ASSUME:  first entry of d an ideal containing numbers and second entry of d 
+         a list of same length containing int vectors with positive entries
+OUTPUT:  a jordan matrix with eigenvalues from the first entry of d 
+         and jordan block sizes from the second entry of d
+EXAMPLE: example jordanmatrix; shows an example
+{
+  def eM,bM=d[1..2];
+
+  int i,j,n;
+  for(i=size(bM);i>=1;i--)
   {
     for(j=size(bM[i]);j>=1;j--)
+    {
+      n=n+bM[i][j];
+    }
+  }
+  matrix M[n][n];
+
+  int k,l;
+  for(i,l,M=1,1,0;i<=size(eM);i++)
+  {
+    for(j=1;j<=size(bM[i]);j++)
     {
       for(k=bM[i][j];k>=2;k,l=k-1,l+1)
@@ -119 +178 @@
   echo=2;
   ring r;
+  list d=ideal(2,3),list(intvec(1,1,2,2,2),intvec(3,4,5));
+  d;
+  print(jordanmatrix(d));
+}
+///////////////////////////////////////////////////////////////////////////////
+
+proc jordanform(matrix M)
+INPUT:   square matrix M with factorizable characteristic polynomial 
+         of the constant part M0
+OUTPUT:  nothing, if the assumptions are not fulfilled, else,
+         a jordan matrix of M0
+EXAMPLE: example jordanform; shows an example
+{
+  return(jordanmatrix(jordandata(M)));
+}
+example
+{
+  "EXAMPLE:";
+  echo=2;
+  LIB "random.lib";
+  ring r;
+  list d=ideal(2,3),list(intvec(1,1,2,2,2),intvec(3,4,5));
+  matrix M=jordanmatrix(d);
+  print(M);
+  int n=nrows(M);
+  matrix U;
+  while(det(U)==0)
+  {
+    U=randommat(n,n,ideal(1)); 
+  }
+  matrix invU=lift(U,freemodule(n));
+  M=invU*M*U;
+  print(M);
+  print(jordanform(M));
+}
+///////////////////////////////////////////////////////////////////////////////
+
+proc jordanbasis(matrix M)
+INPUT:   square matrix M with factorizable characteristic polynomial 
+         of the constant part M0
+OUTPUT:  nothing, if the assumptions are not fulfilled, else,
+         a module containing a basis with respect to which M0 is in
+         jordan normal form
+EXAMPLE: example jordanbasis; shows an example
+{
+  int n=nrows(M);
+  if(n!=ncols(M))
+  {
+    "//no square matrix";
+    return();
+  }
+
+  def br=basering;
+  map zero=br,0;
+  M=zero(M);
+
+  changeord("pr","dp");
+  matrix M=imap(br,M);
+
+  list l=factorize(det(M-var(1)*freemodule(n)),2);
+  def eM,mM=l[1..2];
+
   int i;
-  matrix M[10][10];
-  for(i=1;i<=4;i++) {M[i,i]=2;}
-  for(i=1;i<=3;i++) {M[i,i+1]=1;}
-  for(i=5;i<=10;i++) {M[i,i]=3;}
-  for(i=5;i<=6;i++) {M[i,i+1]=1;}
-  for(i=8;i<=8;i++) {M[i,i+1]=1;}
-  print(M);
-  jordan(M,1);
-  print(jordan(M));
-}
-///////////////////////////////////////////////////////////////////////////////
+  for(i=size(eM);i>=1;i--)
+  {
+    if(deg(eM[i])>1)
+    {
+      kill pr;
+      "//unable to factorize characteristic polynomial";
+      return();
+    }
+  }
+
+  map inv=pr,-var(1);
+  eM=simplify(inv(eM),1);
+  setring br;
+  map zero=pr,0;
+  ideal eM=zero(eM);
+  kill pr;
+
+  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;
+
+  int k,l;
+  matrix Mi,Ni;
+  matrix E=freemodule(n);
+  module V,K1,K2,sNi;
+  matrix v[n][1];
+  list K;
+  for(i=ncols(eM);i>=1;i--)
+  {
+    Mi=M-eM[i]*E;
+    K=list();
+
+    for(Ni,sNi=Mi,0;size(sNi)<mM[i];Ni=Ni*Mi)
+    {
+      sNi=syz(Ni);
+      K=K+list(sNi);
+    }
+    if(size(K)>1)
+    {
+      K1=K[1];
+      K[1]=interred(reduce(K[1],std(module(Mi*K[2]))));
+      for(j=2;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>=1;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;
+        }
+      }
+    }
+  }
+  return(V);
+}
+example
+{
+  "EXAMPLE:";
+  echo=2;
+  LIB "random.lib";
+  ring r;
+  list d=ideal(2,3),list(intvec(1,1,2,2,2),intvec(3,4,5));
+  matrix M=jordanmatrix(d);
+  print(M);
+  int n=nrows(M);
+  matrix U;
+  while(det(U)==0)
+  {
+    U=randommat(n,n,ideal(1)); 
+  }
+  matrix invU=lift(U,freemodule(n));
+  M=invU*M*U;
+  print(M);
+  U=jordanbasis(M);
+  invU=lift(U,freemodule(n));
+  print(invU*M*U);
+}
+///////////////////////////////////////////////////////////////////////////////