source: git/Singular/LIB/jacobson.lib @ 1f190d7

//////////////////////////////////////////////////////////////////////////////
version="$Id$";
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,
@*           Viktor Levandovskyy,      levandov@math.rwth-aachen.de

THEORY: We work over a ring R, that is an Euclidean 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 have two variables, say x and d.
@* We treat then the basering as the Ore localization of R
@* with respect to the mult. closed set S = K[x] without 0.
@* Thus, we treat basering as principal ideal ring with d a polynomial
@* variable and x an invertible one.
@* 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.
@*
@* We call a square matrix D as before 'a weak Jacobson normal form of M'.
@* It is known, that over the first rational Weyl algebra K(x)<d>, D can be further
@* transformed into a diagonal matrix (1,1,...,1,f,0,..,0), where f is in K(x)<d>. We call
@* such a form of D the strong Jacobson normal form. The existence of strong form
@* in not guaranteed if one works with algebra, which is not rational Weyl algebra.


REFERENCES:
@* [1] N. Jacobson, 'The theory of rings', AMS, 1943.
@* [2] 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.
@* [3] V. Levandovskyy, K. Schindelar 'Computing Jacobson normal form using Groebner bases',
@* to appear in Journal of Symbolic Computation, 2010.

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
divideUnits(L);     create ones out of units in the output of smith or jacobson


SEE ALSO: control_lib
";

  LIB "poly.lib";
  LIB "involut.lib"; // involution
  LIB "dmodapp.lib"; // engine
  LIB "qhmoduli.lib"; // Min

proc tstjacobson()
{
  /* tests all procs for consistency */
  example divideUnits;
  example smith;
  example jacobson;
}

proc divideUnits(list L)
"USAGE: divideUnits(L); list L
RETURN: matrix or list of matrices
ASSUME: L is an output of @code{smith} or a @code{jacobson} procedures, that is
@* either L contains one rectangular matrix with elements only on the main diagonal
@* or L consists of three matrices, where L[1] and L[3] are square invertible matrices
@* while L[2] is a rectangular matrix with elements only on the main diagonal
PURPOSE: divide out units from the diagonal and reflect this in transformation matrices
EXAMPLE: example divideUnits; shows examples
"
{
  // check assumptions
  int s = size(L);
  if ( (s!=1) && (s!=3) )
  {
    ERROR("The list has neither 1 nor 3 elements");
  }
  // check sizes of matrices

  if (s==1)
  {
    list LL; LL[1] = L[1]; LL[2] = L[1];
    L = LL;
  }


  // divide units out
  matrix M = L[2];
  int ncM = ncols(M);   int nrM = nrows(M);
  //  matrix TM[nrM][nrM]; // m times m matrix
  matrix TM = matrix(freemodule(nrM));
  int i; int nsize = Min(intvec(nrM,ncM));
  poly p; number n; intvec lexp;

  for(i=1; i<=nsize; i++)
  {
    p = M[i,i];
    lexp = leadexp(p);
    //    TM[i,i] = 1;
    // check: is p a unit?
    if (p!=0)
    {
      if ( lexp == 0)
      {
        // hence p = n*1
        n = leadcoef(p);
        TM[i,i] = 1/n;
      }
    }
  }
  int ppl = printlevel-voice+2;
  dbprint(ppl,"divideUnits: extra transformation matrix is: ");
  dbprint(ppl, TM);
  L[2] = TM*L[2];
  if (s==3)
  {
    L[1] = TM*L[1];
    return(L);
  }
  else
  {
    return(L[2]);
  }
}
example
{ "EXAMPLE:"; echo = 2;
  ring R=(0,m,M,L1,L2,m1,m2,g), D, lp; // two pendula example
  matrix P[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(P,1);  // returns a list with 3 entries
  print(s[2]); // a diagonal form, close to the Smith form
  print(s[1]); // U, left transformation matrix
  list t = divideUnits(s);
  print(t[2]); // the Smith form of the matrix P
  print(t[1]); // U', modified left transformation matrix
}

proc smith(matrix MA, list #)
"USAGE: smith(M[, eng1, eng2]);  M matrix, eng1 and eng2 are optional integers
RETURN: matrix or list of matrices, depending on arguments
ASSUME: Basering is a commutative polynomial ring in one variable
PURPOSE: compute the Smith Normal Form of M with (optionally) transformation matrices
THEORY: Groebner bases are used for the Smith form like in [2] and [3].
NOTE: By default, just the Smith normal form of M is returned.
@* If the optional integer @code{eng1} is non-zero, the list {U,D,V} is returned
@* where U*M*V = D and the diagonal field entries of D are not normalized.
@* The normalization of the latter can be done with the 'divideUnits' procedure.
@* U and V above are square unimodular (invertible) matrices.
@* Note, that the procedure works for a rectangular matrix M.
@*
@* The optional integer @code{eng2} determines the Groebner basis engine:
@* 0 (default) ensures the use of 'slimgb' , otherwise 'std' is used.
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
SEE ALSO: divideUnits, jacobson
"
{
  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;
  if ( size(#)>0 )
  {
    eng=1;
    if (typeof(#[1])=="int")
    {
      eng=int(#[1]); // zero can also happen
    }
  }
  if (size(#)==2)
  {
  BASIS=#[2];
  }
  else{BASIS=0;}


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

  //generate a module consisting of the columns of MA
  module m=MA[1];
  int i;
  for(i=2;i<=ROW;i++)
  {
    m=m,MA[i];
  }

  //if MA eqauls the zero matrix give back MA
  if ( (size(module(m))==0) and (size(transpose(module(m)))==0) )
  {
    module L=freemodule(COL);
    matrix LM=L;
    L=freemodule(ROW);
    matrix RM=L;
    list RUECK=RM;
    RUECK=insert(RUECK,MA);
    RUECK=insert(RUECK,LM);
    return(RUECK);
  }

  if(eng==1)
  {
    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;
  list s=smith(m,1);
  print(s[2]);  // non-normalized Smith form of m
  print(s[1]*m*s[3] - s[2]); // check U*M*V = D
  list t = divideUnits(s);
  print(t[2]); // the Smith form of m
}

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


// VL : engine replaced by the one from dmodapp.lib
// cases are changed

// 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 (non-commutative) ring in two variables.
PURPOSE: compute a weak Jacobson normal form of M over the basering
THEORY: Groebner bases and involutions are used, following [3]
NOTE: A list L of matrices {U,D,V} is returned. That is L[1]*M*L[3]=L[2], 
@*      where L[2] is a diagonal matrix and 
@*      L[1], L[3] are square invertible polynomial (unimodular) matrices.
@*      Note, that M can be rectangular.
@* The optional integer @code{eng2} determines the Groebner basis engine:
@* 0 (default) ensures the use of 'slimgb' , otherwise 'std' is used.
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
SEE ALSO: divideUnits, smith
"
{
  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 )
  {
    B=1;
    if (typeof(#[1])=="int")
    {
    B=int(#[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
  map MAP=R,var(1),var(2);
  matrix M=MAP(MA);

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

  //back to the ring
  setring R;
  map MAPR=r,var(1),var(2);
  module ma=MAPR(m);
  matrix MAA=ma;
  module TL=MAPR(TrafoL);
  module TR=MAPR(TrafoR);
  matrix TRR=TR;
  matrix CON=divideByContent(MAA);

list RUECK=CON*TL, CON*MAA, TRR;
return(RUECK);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,d),Dp;
  def R = nc_algebra(1,1);   setring R; // the 1st Weyl algebra
  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 for m
  print(J[1]*m*J[3] - J[2]); // check that U*M*V = D
  /*   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 terminated /////////////////
      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 consumig
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;

// example from the paper:
    ring w = 0,(x,d),Dp;
    def W=nc_algebra(1,1);
    setring W;
    matrix m[2][2]=d^2-1,d+1,d^2+1,d-x;
    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]);

    ring w2 = 0,(x,s),Dp;
    def W2=nc_algebra(1,s);
    setring W2;
    poly d = s;
    matrix m[2][2]=d^2-1,d+1,d^2+1,d-x;
    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]);
// here, both JNFs are cyclic

// another example from the paper:
    ring w = 0,(x,d),Dp;
    def W=nc_algebra(1,1);
    setring W;
   matrix m[2][2]=-x*d+1, x^2*d, -d, x*d+1;
    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]);

    ring w2 = 0,(x,s),Dp;
    def W2=nc_algebra(1,s);
    setring W2;
    poly d = s;
   matrix m[2][2]=-x*d+1, x^2*d, -d, x*d+1;
    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 example from the paper, also Middeke
   ring w = (0,y),(x,d),Dp;
   def W=nc_algebra(1,1);
   setring W;
   matrix m[2][2]=y^2*d^2+d+1, 1, x*d, x^2*d^2+d+y;
   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]);


*/