Ticket #57: jacobson.lib.diff.txt

//////////////////////////////////////////////////////////////////////////////
version="$Id: jacobson.lib,v 1.9 2009/02/12 20:27:17 levandov Exp $";
category="System and Control Theory";
info="
LIBRARY: jacobson.lib     Algorithms for Smith and Jacobson Normal Form
AUTHOR: Kristina Schindelar,     Kristina.Schindelar@math.rwth-aachen.de

THEORY: We work over a ring R, that is a principal ideal domain.
@*   If R is commutative, we suppose R to be a polynomial ring in one variable.
@* If R is non-commutative, we suppose R to be in two variables, say x and d.
@* We treat then the basering as principal ideal ring with d a polynomial
@* variable and x an invertible one. That is, we work in the Ore localization of R
@* with respect to the mult. closed set S = K[x] without 0.
@* Note, that in computations no division by x will actually happen.
@*   Given a rectangular matrix M over R, one can compute unimodular (that is invertible)
@* square matrices U and V, such that U*M*V=D is a diagonal matrix.
@*   Depending on the ring, the diagonal entries of D have certain properties.

REFERENCES:
@* N. Jacobson, 'The theory of rings', AMS, 1943.
@* Manuel Avelino Insua Hermo, 'Varias perspectives sobre las bases de Groebner: Forma normal de Smith, Algorithme de Berlekamp y algebras de Leibniz'. PhD thesis, Universidad de Santiago de Compostela, 2005.


PROCEDURES:
smith(M[,eng1,eng2]);  compute the Smith Normal Form of M over commutative ring
jacobson(M[,eng]);       compute a weak Jacobson Normal Form of M over non-commutative ring

SEE ALSO: control_lib
";

  LIB "poly.lib";
  LIB "involut.lib";

proc smith(matrix MA, list #)
"USAGE: smith(M[, eng1, eng2]);  M matrix, eng1 and eng2 are optional integers
RETURN: matrix or list
ASSUME: The current ring is assumed to be the commutative polynomial ring in
        one variable
PURPOSE: compute the Smith Normal Form of M with transformation matrices (optionally)
NOTE: If the optional integer eng1 is non-zero, returns the list {U,D,V},
where U*M*V = D and the diagonal field entries of D are not normalized.
@*    Otherwise only the matrix D, that is the Smith Normal Form of M, is returned.
@*    Here, U and V are square unimodular (invertible) matrices. The procedure works for rectangular matrix M.
@* The optional integer eng2 determines the engine, that computes the Groebner basis:
@* 0 means 'std' (default), 1 means 'groebner' and 2 means 'slimgb'.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@* if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example smith; shows examples
"
{
  def R = basering;
  // check assume: R must be a polynomial ring in 1 variable
  setring R;
  if (nvars(R) >1 )
  {
    ERROR("Basering must have exactly one variable");
  }

  int eng = 0;
  int BASIS = 0;

  if ( size(#)>0 )
  {
    if (typeof(#[1])=="int")
    {
      eng=#[1]; // zero can also happen
  }

    if (size(#) > 1)
    {
      if (typeof(#[2])=="int")
  {
  BASIS=#[2];
  }
    } 
  }


  int ROW=ncols(MA);
  int COL=nrows(MA);

  //generate a module consisting of the columns of MA
  module m = MA;
  int i;

  //if MA eqauls the zero matrix give back MA
  if ( (size(m)==0) and (size(transpose(m))==0) )
  {
    if( eng > 0 )
    {
      return(list( matrix(freemodule(COL)), MA, matrix(freemodule(ROW)) ));
    }
   
    return(MA);
  }

  if(eng > 0)
  {
    list rueckLI=diagonal_with_trafo(R,MA,BASIS);
    list rueckLII=divisibility(rueckLI[2]);
    matrix CON=divideByContent(rueckLII[2]);
    list rueckL=CON*rueckLII[1]*rueckLI[1], CON*rueckLII[2], rueckLI[3]*rueckLII[3];
    return(rueckL);
  }
  else
  {
    matrix rueckm=diagonal_without_trafo(R,MA,BASIS);
    list rueckL=divisibility(rueckm);
    matrix CON=divideByContent(rueckm);
    rueckm=CON*rueckL[2];
    return(rueckm);
  }
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = 0,x,Dp;
  matrix m[3][2]=x, x^4+x^2+21, x^4+x^2+x, x^3+x, 4*x^2+x, x;
  print(m); // M
  print(smith(m)); // Smith Normal Form of M
  list s=smith(m,1);
  // TODO: BUG THE FOLLOWING GIVES A DIFFERENT RESULT TO ABOVE!
  print(S[2]); // Smith Normal Form D of M 
  // TODO: ARE U AND V REALLY UNIMODULAR?
  S[1]; // U
  det(S[1]);
  S[3]; // V
  det(S[3]);
  print(S[1]*m*S[3] - S[2]); // check that U*M*V = D
  print(smith(matrix(0))); // Smith Normal Form of 0

}

static proc diagonal_with_trafo( R, matrix MA, int B)
{

int ppl = printlevel-voice+2;

  int BASIS=B;
  int ROW=ncols(MA);
  int COL=nrows(MA);
  module m=MA[1];
  int i,j,k;
  for(i=2;i<=ROW;i++)
  {
    m=m,MA[i];
  }


  //add zero rows or columns
  //add zero rows or columns
  int adrow=0;
  for(i=1;i<=COL;i++)
  {
  k=0;
  for(j=1;j<=ROW;j++)
  {
  if(MA[i,j]!=0){k=1;}
  }
  if(k==0){adrow++;}
  }

    m=transpose(m);
    for(i=1;i<=adrow;i++){m=m,0;}
    m=transpose(m);

  list RINGLIST=ringlist(R);
  list o="C",0;
  list oo="lp",1;
  list ORD=o,oo;
  RINGLIST[3]=ORD;
  def r=ring(RINGLIST);
  setring r;
  //fix the required ordering
  map MAP=R,var(1);
  module m=MAP(m);

  int flag=1;  ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0)

    module TrafoL=freemodule(COL);
    module TrafoR=freemodule(ROW);
    module EXL=freemodule(COL); //because we start with transpose(m)
    module EXR=freemodule(ROW);

    option(redSB);
    option(redTail);
    module STD_EX;
    module Trafo;

    int s,st,p,ff;

    module LT,TSTD;
    module STD=transpose(m);
    int finish=0;
    int fehlpos;
    list pos;
    matrix END;
    matrix ENDSTD;
    matrix STDFIN;
    STDFIN=STD;
    list COMPARE=STDFIN;

    while(finish==0)
    {
      dbprint(ppl,"Going into the while cycle");

      if(flag mod 2==1)
      {
        STD_EX=EXL,transpose(STD);
      }
      else
      {
        STD_EX=EXR,transpose(STD);
      }
        dbprint(ppl,"Computing Groebner basis: start");
        dbprint(ppl-1,STD);
      STD=engine(STD,BASIS);
        dbprint(ppl,"Computing Groebner basis: finished");
        dbprint(ppl-1,STD);

      if(flag mod 2==1){s=ROW; st=COL;}else{s=COL; st=ROW;}
        dbprint(ppl,"Computing Groebner basis for transformation matrix: start");
        dbprint(ppl-1,STD_EX);
      STD_EX=transpose(STD_EX);
      STD_EX=engine(STD_EX,BASIS);
        dbprint(ppl,"Computing Groebner basis for transformation matrix: finished");
        dbprint(ppl-1,STD_EX);

      //////// split STD_EX in STD and the transformation matrix
      STD_EX=transpose(STD_EX);
      STD=STD_EX[st+1];
      LT=STD_EX[1];

      ENDSTD=STD_EX;
      for(i=2; i<=ncols(ENDSTD); i++)
      {
        if (i<=st)
        {
         LT=LT,STD_EX[i];
        }
        if (i>st+1)
        {
          STD=STD,STD_EX[i];
        }
      }

      STD=transpose(STD);
      LT=transpose(LT);

      ////////////////////// compute the transformation matrices

      if (flag mod 2 ==1)
      {
        TrafoL=transpose(LT)*TrafoL;
      }
      else
      {
        TrafoR=TrafoR*LT;
      }


      STDFIN=STD;
      /////////////////////////////////// check if the D-th row is finished ///////////////////////////////////
      COMPARE=insert(COMPARE,STDFIN);
             if(size(COMPARE)>=3)
       {
        if(STD==COMPARE[3]){finish=1;}
       }
////////////////////////////////// change to the opposite module
       TSTD=transpose(STD);
       STD=TSTD;
       flag++;
      dbprint(ppl,"Finished one while cycle");
    }


   if (flag mod 2!=0) { STD=transpose(STD); }


   //zero colums to the end
    matrix STDMM=STD;
    pos=list();
    fehlpos=0;
    while( size(STD)+fehlpos-ncols(STDMM) < 0)
    {
      for(i=1; i<=ncols(STDMM); i++)
      {
        ff=0;
        for(j=1; j<=nrows(STDMM); j++)
        {
          if (STD[j,i]==0) { ff++; }
        }
        if(ff==nrows(STDMM))
        {
          pos=insert(pos,i); fehlpos++;
        }
      }
    }
    int fehlposc=fehlpos;
    module SORT;
    for(i=1; i<=fehlpos; i++)
    {
      SORT=gen(2);
      for (j=3;j<=ROW;j++)
      {
        SORT=SORT,gen(j);
      }
      SORT=SORT,gen(1);
      STD=STD*SORT;
      TrafoR=TrafoR*SORT;
    }

    //zero rows to the end
    STDMM=transpose(STD);
    pos=list();
    fehlpos=0;
    while( size(transpose(STD))+fehlpos-ncols(STDMM) < 0)
    {
      for(i=1; i<=ncols(STDMM); i++)
      {
        ff=0; for(j=1; j<=nrows(STDMM); j++)
        {
           if(transpose(STD)[j,i]==0){ ff++;}}
           if(ff==nrows(STDMM)) { pos=insert(pos,i); fehlpos++; }
      }
    }
    int fehlposr=fehlpos;

    for(i=1; i<=fehlpos; i++)
    {
      SORT=gen(2);
      for(j=3;j<=COL;j++){SORT=SORT,gen(j);}
      SORT=SORT,gen(1);
      SORT=transpose(SORT);
      STD=SORT*STD;
      TrafoL=SORT*TrafoL;
    }

    setring R;
    map MAPinv=r,var(1);
    module STD=MAPinv(STD);
    module TrafoL=MAPinv(TrafoL);
    matrix TrafoLM=TrafoL;
    module TrafoR=MAPinv(TrafoR);
    matrix TrafoRM=TrafoR;
    matrix STDM=STD;

    //Test
    if(TrafoLM*m*TrafoRM!=STDM){ return(0); }

    list RUECK=TrafoRM;
    RUECK=insert(RUECK,STDM);
    RUECK=insert(RUECK,TrafoLM);
    return(RUECK);
}


static proc divisibility(matrix M)
   {
    matrix STDM=M;
    int i,j;
    int ROW=nrows(M);
    int COL=ncols(M);
    module TrafoR=freemodule(COL);
    module TrafoL=freemodule(ROW);
    module SORT;
    matrix TrafoRM=TrafoR;
    matrix TrafoLM=TrafoL;
    list posdeg0;
    int posdeg=0;
    int act;
    int Sort=ROW;
    if(size(std(STDM))!=0)
    { while( size(transpose(STDM)[Sort])==0 ){Sort--;}}

    for(i=1;i<=Sort ;i++)
    {
      if(leadexp(STDM[i,i])==0){posdeg0=insert(posdeg0,i);posdeg++;}
    }
    //entries of degree 0 at the beginning
    for(i=1; i<=posdeg; i++)
    {
      act=posdeg0[i];
      SORT=gen(act);
      for(j=1; j<=COL; j++){if(j!=act){SORT=SORT,gen(j);}}
      STDM=STDM*SORT;
      TrafoRM=TrafoRM*SORT;
      SORT=gen(act);
      for(j=1; j<=ROW; j++){if(j!=act){SORT=SORT,gen(j);}}
      STDM=transpose(SORT)*STDM;
      TrafoLM=transpose(SORT)*TrafoLM;
    }

    poly G;
    module UNITL=freemodule(ROW);
    matrix GCDL=UNITL;
    module UNITR=freemodule(COL);
    matrix GCDR=UNITR;
    for(i=posdeg+1; i<=Sort; i++)
    {
      for(j=i+1; j<=Sort; j++)
      {
        GCDL=UNITL;
        GCDR=UNITR;
        G=gcd(STDM[i,i],STDM[j,j]);
        ideal Z=STDM[i,i],STDM[j,j];
        matrix T=lift(Z,G);
        GCDL[i,i]=T[1,1];
        GCDL[i,j]=T[2,1];
        GCDL[j,i]=-STDM[j,j]/G;
        GCDL[j,j]=STDM[i,i]/G;
        GCDR[i,j]=T[2,1]*STDM[j,j]/G;
        GCDR[j,j]=T[2,1]*STDM[j,j]/G-1;
        GCDR[j,i]=1;
        STDM=GCDL*STDM*GCDR;
        TrafoLM=GCDL*TrafoLM;
        TrafoRM=TrafoRM*GCDR;
      }
    }
    list RUECK=TrafoRM;
    RUECK=insert(RUECK,STDM);
    RUECK=insert(RUECK,TrafoLM);
    return(RUECK);
}

static proc diagonal_without_trafo( R, matrix MA, int B)
{
  int ppl = printlevel-voice+2;

  int BASIS=B;
  int ROW=ncols(MA);
  int COL=nrows(MA);
  module m=MA[1];
  int i;
  for(i=2;i<=ROW;i++)
  {m=m,MA[i];}


  list RINGLIST=ringlist(R);
  list o="C",0;
  list oo="lp",1;
  list ORD=o,oo;
  RINGLIST[3]=ORD;
  def r=ring(RINGLIST);
  setring r;
  //RICHTIGE ORDNUNG MACHEN
  map MAP=R,var(1);
  module m=MAP(m);

  int flag=1; ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0)


    int act, j, ff;
    option(redSB);
    option(redTail);


    module STD=transpose(m);
    module TSTD;
    int finish=0;
    matrix STDFIN;
    STDFIN=STD;
    list COMPARE=STDFIN;

    while(finish==0)
    {
      dbprint(ppl,"Going into the while cycle");
      dbprint(ppl,"Computing Groebner basis: start");
      dbprint(ppl-1,STD);
      STD=engine(STD,BASIS);
      dbprint(ppl,"Computing Groebner basis: finished");
      dbprint(ppl-1,STD);
      STDFIN=STD;
      /////////////////////////////////// check if the D-th row is finished ///////////////////////////////////
      COMPARE=insert(COMPARE,STDFIN);
             if(size(COMPARE)>=3)
       {
        if(STD==COMPARE[3]){finish=1;}
       }
        ////////////////////////////////// change to the opposite module

        TSTD=transpose(STD);
        STD=TSTD;
        flag++;
        dbprint(ppl,"Finished one while cycle");
    }

    matrix STDMM=STD;
    list pos=list();
    int fehlpos=0;
    while( size(STD)+fehlpos-ncols(STDMM) < 0)
    {
      for(i=1; i<=ncols(STDMM); i++)
      {
        ff=0;
        for(j=1; j<=nrows(STDMM); j++)
        {
          if (STD[j,i]==0) { ff++; }
        }
        if(ff==nrows(STDMM))
        {
          pos=insert(pos,i); fehlpos++;
        }
      }
    }
   int fehlposc=fehlpos;
    module SORT;
    for(i=1; i<=fehlpos; i++)
    {
      SORT=gen(2);
      for (j=3;j<=ROW;j++)
      {
        SORT=SORT,gen(j);
      }
      SORT=SORT,gen(1);
      STD=STD*SORT;
    }

    //zero rows to the end
    STDMM=transpose(STD);
    pos=list();
    fehlpos=0;
    while( size(transpose(STD))+fehlpos-ncols(STDMM) < 0)
    {
      for(i=1; i<=ncols(STDMM); i++)
      {
        ff=0; for(j=1; j<=nrows(STDMM); j++)
        {
           if(transpose(STD)[j,i]==0){ ff++;}}
           if(ff==nrows(STDMM)) { pos=insert(pos,i); fehlpos++; }
      }
    }
    int fehlposr=fehlpos;

    for(i=1; i<=fehlpos; i++)
    {
      SORT=gen(2);
      for(j=3;j<=COL;j++){SORT=SORT,gen(j);}
      SORT=SORT,gen(1);
      SORT=transpose(SORT);
      STD=SORT*STD;
    }

   //add zero rows or columns

    int adrow=COL-size(transpose(STD));
    int adcol=ROW-size(STD);

    for(i=1;i<=adcol;i++){STD=STD,0;}
    STD=transpose(STD);
    for(i=1;i<=adrow;i++){STD=STD,0;}
    STD=transpose(STD);

    setring R;
    map MAPinv=r,var(1);
    module STD=MAPinv(STD);
    matrix STDM=STD;
    return(STDM);
}



static proc engine(module I, int i)
{
  module J;
  if (i==0)
  {
    J = std(I);
  }
  if (i==1)
  {
    J = groebner(I);
  }
  if (i==2)
  {
    J = slimgb(I);
  }
  return(J);
}

proc jacobson(matrix MA, list #)
 "USAGE:  jacobson(M[, eng]);  M matrix, eng an optional int
RETURN: list
ASSUME: Basering is a noncommutative ring in two variables.
PURPOSE: compute a weak Jacobson Normal Form of M over a noncommutative ring
NOTE: A list L of matrices {U,D,V} is returned. That is U*M*V=D, where
@*      D is a diagonal matrix and U, V square invertible (unimodular) matrices.
@*      Note, that M can be rectangular.
@* The optional integer @code{eng} determines the engine, that computes the Groebner basis:
@* 0 means 'std' (default), 1 means 'groebner' and 2 means 'slimgb'.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@* if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example jacobson; shows examples
"
{
  def R = basering;
  // check assume: R must be a polynomial ring in 2 variables
  setring R;
  if ( (nvars(R) !=2 ) )
  {
    ERROR("Basering must have exactly two variables");
  }

  // check if MA is zero matrix and return corr. U,V
  if ( (size(module(MA))==0) and (size(transpose(module(MA)))==0) )
  {
    int nr = nrows(MA);
    int nc = ncols(MA);
    matrix U = matrix(freemodule(nr));
    matrix V = matrix(freemodule(nc));
    list L; L[1]=U; L[2] = MA; L[3] = V;
    return(L);
  }

  int B=0;
  if ( size(#)>0 )
  {
    if (typeof(#[1])=="int")
    {
      B = #[1]; // zero can also happen
    }
  }

  //change ring
  list RINGLIST=ringlist(R);
  list o="C",0;
  intvec v=0,1;
  list oo="a",v;
  v=1,1;
  list ooo="lp",v;
  list ORD=o,oo,ooo;
  RINGLIST[3]=ORD;
  def r=ring(RINGLIST);
  setring r;

  //fix the required ordering
  matrix M=imap(R, MA);

  list T = triangle(M,B);
  module TrafoL = T[1];
  module m      = T[2];
  module TrafoR = T[3];

  //back to the ring
  setring R;

  matrix MAA = imap(r, m);
  matrix CON=divideByContent(MAA);
  matrix TL = imap(r, TrafoL);
  matrix TR = imap(r, TrafoR);
 
  // TODO: CHECK WHETHER HERE SHOULD BE CON*TR!?
  return(list(CON*TL, CON*MAA, TR)); 

}
example
{ "EXAMPLE:"; echo = 2;
  ring r = 0,(x,d),Dp;
  def R=nc_algebra(1,1); setring R;
  R; // the Weyl algebra in x and d
  jacobson(matrix(0));
  
  matrix m[2][2]=d,x,0,d; 
  print(m);
  list J = jacobson(m); // returns a list with 3 entries
  print(J[2]); // a Jacobson Form D
  J[1]; // U
  J[3]; // V

  print(J[1]*m*J[3] - J[2]); // check that U*M*V = D

  matrix m[3][2]=x, x^4+x^2+21, x^4+x^2+x, x^3+x, 4*x^2+x, x;
  print(m); // M
  list JJ = jacobson(m); // returns a list with 3 entries

  // TODO: BUG: WHY THE FOLLOWING JACOBSON FORM IS NOT DIAGONAL???!
  print(JJ[2]); // a Jacobson Form D
  print(JJ[1]*m*JJ[3] - JJ[2]); // check that U*M*V = D
  JJ[1]; // U
  JJ[3]; // V

  list S=smith(m,1);

  print(S[2]); // Smith Normal Form S of M
  print(S[1]*m*S[3] - S[2]); // check that U*M*V = S

  S[1]; // U
  S[3]; // V

  /*   now, let us do the same for the shift algebra  */
  ring r2 = 0,(x,s),Dp;
  def R2 = nc_algebra(1,s);   setring R2; // the 1st shift algebra
  matrix m[2][2] = s,x,0,s; print(m); // matrix of the same for as above
  list J = jacobson(m);
  print(J[2]); // a Jacobson Form D, quite different from above
  print(J[1]*m*J[3] - J[2]); // check that U*M*V = D

}

static proc triangle( matrix MA, int B)
{
 int ppl = printlevel-voice+2;

  map inv=ncdetection();
  int ROW=ncols(MA);
  int COL=nrows(MA);

  //generate a module consisting of the columns of MA
  module m=MA[1];
  int i,j,s,st,p,k,ff,ex, nz, ii,nextp;
  for(i=2;i<=ROW;i++)
  {
    m=m,MA[i];
  }
  int BASIS=B;

  //add zero rows or columns
  int adrow=0;
  for(i=1;i<=COL;i++)
  {
  k=0;
  for(j=1;j<=ROW;j++)
  {
  if(MA[i,j]!=0){k=1;}
  }
  if(k==0){adrow++;}
  }

    m=transpose(m);
    for(i=1;i<=adrow;i++){m=m,0;}
    m=transpose(m);


    int flag=1;  ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0)

    module TrafoL=freemodule(COL);
    module TrafoR=freemodule(ROW);
    module EXL=freemodule(COL); //because we start with transpose(m)
    module EXR=freemodule(ROW);

    option(redSB);
    option(redTail);
    module STD_EX,LT,TSTD, L, Trafo;



    module STD=transpose(m);
    int finish=0;
    list pos, COM, COM_EX;
    matrix END, ENDSTD, STDFIN;
    STDFIN=STD;
    list COMPARE=STDFIN;


  while(finish==0)
    {
      dbprint(ppl,"Going into the while cycle");
      if(flag mod 2==1){STD_EX=EXL,transpose(STD); ex=2*COL;} else {STD_EX=EXR,transpose(STD); ex=2*ROW;}

      dbprint(ppl,"Computing Groebner basis: start");
      dbprint(ppl-1,STD);
      STD=engine(STD,BASIS);
      dbprint(ppl,"Computing Groebner basis: finished");
      dbprint(ppl-1,STD);
      if(flag mod 2==1){s=ROW; st=COL;}else{s=COL; st=ROW;}

      STD_EX=transpose(STD_EX);
      dbprint(ppl,"Computing Groebner basis for transformation matrix: start");
      dbprint(ppl-1,STD_EX);
      STD_EX=engine(STD_EX,BASIS);
      dbprint(ppl,"Computing Groebner basis for transformation matrix: finished");
      dbprint(ppl-1,STD_EX);

      COM=1;
      COM_EX=1;
      for(i=2; i<=size(STD); i++)
         { COM=COM[1..size(COM)],i; COM_EX=COM_EX[1..size(COM_EX)],i; }
      nz=size(STD_EX)-size(STD);

      //zero rows in the begining

       if(size(STD)!=size(STD_EX) )
       {
         for(i=1; i<=size(STD_EX)-size(STD); i++)
         {
           COM_EX=0,COM_EX[1..size(COM_EX)];
         }
       }




       for(i=nz+1; i<=size(STD_EX); i++ )
        {if( leadcoef(STD[i-nz])!=leadcoef(STD_EX[i]) )          {STD[i-nz]=leadcoef(STD_EX[i])*STD[i-nz];}
        }

        //assign the zero rows in STD_EX

       for (j=2; j<=nz; j++)
       {
        p=0;
        i=1;
        while(STD_EX[j-1][i]==0){i++;};
        p=i-1;
        nextp=1;
        k=0;
        i=1;
        while(STD_EX[j][i]==0 and i<=p)
        { k++; i++;}
        if (k==p){ COM_EX[j]=-1; }
       }

      //assign the zero rows in STD
      for (j=2; j<=size(STD); j++)
       {
        i=size(transpose(STD));
        while(STD[j-1][i]==0){i--;}
        p=i;
        i=size(transpose(STD[j]));
        while(STD[j][i]==0){i--;}
        if (i==p){ COM[j]=-1; }
       }

       for(j=1; j<=size(COM); j++)
        {
        if(COM[j]<0){COM=delete(COM,j);}
        }

      for(i=1; i<=size(COM_EX); i++)
        {ff=0;
         if(COM_EX[i]==0){ff=1;}
         else
          {  for(j=1; j<=size(COM); j++)
               { if(COM_EX[i]==COM[j]){ff=1;} }
          }
         if(ff==0){COM_EX[i]=-1;}
        }

      L=STD_EX[1];
      for(i=2; i<=size(COM_EX); i++)
       {
        if(COM_EX[i]!=-1){L=L,STD_EX[i];}
       }

      //////// split STD_EX in STD and the transformation matrix

       L=transpose(L);
       STD=L[st+1];
       LT=L[1];


       for(i=2; i<=s+st; i++)
        {
         if (i<=st)
          {
           LT=LT,L[i];
          }
         if (i>st+1)
          {
          STD=STD,L[i];
          }
        }

       STD=transpose(STD);
       STDFIN=matrix(STD);
       COMPARE=insert(COMPARE,STDFIN);
       LT=transpose(LT);

      ////////////////////// compute the transformation matrices

       if (flag mod 2 ==1)
        {
         TrafoL=transpose(LT)*TrafoL;
        }
       else
        {
         TrafoR=TrafoR*involution(LT,inv);
        }


      ///////////////////////// check whether the alg termined /////////////////
      if(size(COMPARE)>=3)
       {
        if(STD==COMPARE[3]){finish=1;}
       }
       ////////////////////////////////// change to the opposite module
       TSTD=transpose(STD);
       STD=involution(TSTD,inv);
       flag++;
       dbprint(ppl,"Finished one while cycle");
    }

if (flag mod 2 ==0){ STD = involution(STD,inv);} else { STD = transpose(STD);  }

    list REVERSE=TrafoL,STD,TrafoR;
    return(REVERSE);
}

static proc divideByContent(matrix M)
{
//find first entrie not equal to zero
int i,k;
k=1;
vector CON;
   for(i=1;i<=ncols(M);i++)
   {
    if(leadcoef(M[i])!=0){CON=CON+leadcoef(M[i])*gen(k); k++;}
   }
poly con=content(CON);
matrix TL=1/con*freemodule(nrows(M));
return(TL);
}


/////interesting examples for smith////////////////

/*

//static proc Ex_One_wheeled_bicycle()
{
ring RA=(0,m), D, lp;
matrix bicycle[2][3]=(1+m)*D^2, D^2, 1, D^2, D^2-1, 0;
list s=smith(bicycle, 1,0);
print(s[2]);
print(s[1]*bicycle*s[3]-s[2]);
}


//static proc Ex_RLC-circuit()
{
ring RA=(0,m,R1,R2,L,C), D, lp;
matrix  circuit[2][3]=D+1/(R1*C), 0, -1/(R1*C), 0, D+R2/L, -1/L;
list s=smith(circuit, 1,0);
print(s[2]);
print(s[1]*circuit*s[3]-s[2]);
}


//static proc Ex_two_pendula()
{
ring RA=(0,m,M,L1,L2,m1,m2,g), D, lp;
matrix pendula[3][4]=m1*L1*D^2,m2*L2*D^2,(M+m1+m2)*D^2,-1,m1*L1^2*D^2-m1*L1*g,0,m1*L1*D^2,0,0,
m2*L2^2*D^2-m2*L2*g,m2*L2*D^2,0;
list s=smith(pendula, 1,0);
print(s[2]);
print(s[1]*pendula*s[3]-s[2]);
}

//static proc Ex_linerized_satellite_in_a_circular_equatorial_orbit()
{
ring RA=(0,m,omega,r,a,b), D, lp;
matrix satellite[4][6]=
D,-1,0,0,0,0,
-3*omega^2,D,0,-2*omega*r,-a/m,0,
0,0,D,-1,0,0,
0,2*omega/r,0,D,0,-b/(m*r);
list s=smith(satellite, 1,0);
print(s[2]);
print(s[1]*satellite*s[3]-s[2]);
}

//static proc Ex_flexible_one_link_robot()
{
ring RA=(0,M11,M12,M13,M21,M22,M31,M33,K1,K2), D, lp;
matrix robot[3][4]=M11*D^2,M12*D^2,M13*D^2,-1,M21*D^2,M22*D^2+K1,0,0,M31*D^2,0,M33*D^2+K2,0;
list s=smith(robot, 1,0);
print(s[2]);
print(s[1]*robot*s[3]-s[2]);
}

*/


/////interesting examples for jacobson////////////////

/*

//static proc Ex_compare_output_with_maple_package_JanetOre()
{
    ring w = 0,(x,d),Dp;
    def W=nc_algebra(1,1);
    setring W;
    matrix m[3][3]=[d2,d+1,0],[d+1,0,d3-x2*d],[2d+1, d3+d2, d2];
    list J=jacobson(m,0);
    print(J[1]*m*J[3]);
    print(J[2]);
    print(J[1]);
    print(J[3]);
    print(J[1]*m*J[3]-J[2]);
}

// modification for shift algebra
{
    ring w = 0,(x,s),Dp;
    def W=nc_algebra(1,s);
    setring W;
    matrix m[3][3]=[s^2,s+1,0],[s+1,0,s^3-x^2*s],[2*s+1, s^3+s^2, s^2];
    list J=jacobson(m,0);
    print(J[1]*m*J[3]);
    print(J[2]);
    print(J[1]);
    print(J[3]);
    print(J[1]*m*J[3]-J[2]);
}

//static proc Ex_cyclic()
{
    ring w = 0,(x,d),Dp;
    def W=nc_algebra(1,1);
    setring W;
    matrix m[3][3]=d,0,0,x*d+1,d,0,0,x*d,d;
    list J=jacobson(m,0);
    print(J[1]*m*J[3]);
    print(J[2]);
    print(J[1]);
    print(J[3]);
    print(J[1]*m*J[3]-J[2]);
}

// modification for shift algebra: gives a module with ann = s^2
{
    ring w = 0,(x,s),Dp;
    def W=nc_algebra(1,s);
    setring W;
    matrix m[3][3]=s,0,0,x*s+1,s,0,0,x*s,s;
    list J=jacobson(m,0);
    print(J[1]*m*J[3]);
    print(J[2]);
    print(J[1]);
    print(J[3]);
    print(J[1]*m*J[3]-J[2]);
}

// yet another modification for shift algebra: gives a module with ann = s
// xs+1 is replaced by x*s+s
{
    ring w = 0,(x,s),Dp;
    def W=nc_algebra(1,s);
    setring W;
    matrix m[3][3]=s,0,0,x*s+s,s,0,0,x*s,s;
    list J=jacobson(m,0);
    print(J[1]*m*J[3]);
    print(J[2]);
    print(J[1]);
    print(J[3]);
    print(J[1]*m*J[3]-J[2]);
}

// variations for above setup with shift algebra:

// easy but instructive one: see the difference to Weyl case!
matrix m[2][2]=s,x,0,s; print(m);

// more interesting:
matrix n[3][3]=s,x,0,0,s,x,s,0,x;

// things blow up
matrix m[3][4] = s,x^2*s,x^3*s,s*x^2,s*x+1,(x+1)^3, (x+s)^2, x*s;

// this one is quite nasty and time consuming
matrix m[3][4] = s,x^2*s,x^3*s,s*x^2,s*x+1,(x+1)^3, (x+s)^2, x*s,x,x^2,x^3,s;

*/