Ticket #57: jacobson.lib

//////////////////////////////////////////////////////////////////////////////
version="$Id: jacobson.lib,v 1.4 2008/12/01 20:51:16 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 principal ideal domain R. We suppose R to be either 
@*  the commutative polynomial ring in one variable or the Ore localization of 
@*  the first Weyl algebra with respect to S = K[x]\{0} (i.e. K(x)<D|D*x=x*D+1>)
@*  Under these assumptions on R, given a rectangular matrix M over R, one can 
@*  compute unimodular matrices U, V such that U*M*V=D is a diagonal matrix. 
@*  Depending on the ring, the diagonal entries of D have certain properties.
@*  Over a commutative ring in one variable, the matrix D is unique and called 
@*  the Smith Normal Form of M, whereas in the non-commutative case D is not 
@*  unique and a weak Jacobson Normal Form of M is computed.

TODO: WHAT ABOUT UNIQUENESS OF JACOBSON NORMAL FORM? WHICH ONE IS RETURNED?

PROCEDURES:
smith(M[,eng1,eng2]);  computes the Smith Normal Form
jacobson(M[,eng]);     computes a weak Jacobson Normal Form

SEE ALSO: control_lib
";

  LIB "poly.lib"; // TODO: WHAT FOR?
  LIB "involut.lib"; // TODO: WHAT FOR?

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

// TODO: IT IS BETTER TO WORK WITH MODULES INSTEAD OF MATRICES!!!
// PLEASE SWITCH MATRICES TO MODULES!

// TODO: IS THE ALGORITHM WELL KNOWN OR IS IT DUE TO SOME BOOK??? 
// ADD REFERENCE IS SO.
proc smith(matrix MA, list #)
"USAGE: smith(M[, eng1, eng2]);  M matrix, eng1 and eng2 are optional integers
RETURN: matrix or list, either a matrix D or a list containing matrices {U,D,V},
       such that U*M*V = D, where D is diagonal (the Smith Normal Form of M) 
       and the matrices U and V are unimodular.
ASSUME: The current ring is assumed to be the commutative polynomial ring in 
        one variable
NOTE: If the optional integer eng1 is non-zero, returns the list {U,D,V}, 
@*    otherwise only the matrix, the Smith Normal Form of M is returned.
@*    The optional integer eng2 determines the engine, that computes the GB:
@*    0 means 'std' (default), 1 means 'groebner' and 2 means 'slimgb'.
@* 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;

  // TODO: ADD ASSUMPTION CHECKS!!!

  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) )
  {
    return(list( matrix(freemodule(COL)), MA, matrix(freemodule(ROW)) ));
  }

  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

}

// TODO: THIS PROCEDURE SEEMS TOO COMPLICATED, THUS PLEASE
// TODO: 1. ADD HELP WITH ASSUMPTIONS AND RETURN DESCRIPTION
// TODO: 2. ADD STEP-BY-STEP COMMENTS INTO THE CODE

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);
}



// TODO: THIS PROCEDURE SEEMS TOO COMPLICATED, THUS PLEASE
// TODO: 1. ADD HELP WITH ASSUMPTIONS AND RETURN DESCRIPTION
// TODO: 2. ADD STEP-BY-STEP COMMENTS INTO THE CODE

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);
}


// TODO: THIS PROCEDURE SEEMS TOO COMPLICATED, THUS PLEASE
// TODO: 1. ADD HELP WITH ASSUMPTIONS AND RETURN DESCRIPTION
// TODO: 2. ADD STEP-BY-STEP COMMENTS INTO THE CODE

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);
}

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


// TODO: IS THE ALGORITHM WELL KNOWN OR IS IT DUE TO SOME BOOK??? 
// ADD REFERENCE IF SO.
// TODO: ADD A NOTE ABOUT THE CHOSEN JACOBSON NORMAL FORM!

proc jacobson(matrix MA, list #)
"USAGE:  jacobson(M[, eng]);  M matrix, eng an optional integer
RETURN: list, a list containing matrices {U,D,V}, such that U*M*V = D, 
        where U and V are unimodular, D is diagonal.
ASSUME: The current ring is assumed to be the Ore localization of the first 
        Weyl algebra with respect to S = K[x]\{0} (i.e. K(x)<D|D*x=x*D+1>).
NOTE: The optional integer eng2 determines the engine, that computes the GB:
@*    0 means 'std' (default), 1 means 'groebner' and 2 means 'slimgb'.
@* 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;

  // TODO: ADD ASSUMPTION CHECKS!!!

  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=1,0;
  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)); // TODO: BUG HERE!!!

  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
  // TODO: ARE U AND V REALLY UNIMODULAR?
  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

  // TODO: ARE U AND V REALLY UNIMODULAR?
  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

  // TODO: ARE U AND V REALLY UNIMODULAR?
  S[1]; // U
  S[3]; // V
}


// TODO: THIS PROCEDURE SEEMS TOO COMPLICATED, THUS PLEASE
// TODO: 1. ADD HELP WITH ASSUMPTIONS AND RETURN DESCRIPTION
// TODO: 2. ADD STEP-BY-STEP COMMENTS INTO THE CODE

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(RA,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(RA,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(RA,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(RA,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(RA,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(W,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(W,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]);
}