source: git/Singular/LIB/multigrading.lib @ 087946

///////////////////////////////////////////////////////////////////
Version="$Id: multigrading.lib 0.1 2010-02-13$";
category="Combinatorial Commutative Algebra";
info="
LIBRARY:  multigrading.lib          Multigraded Rings

AUTHORS:  Rene Birkner, rbirkner@math.fu-berlin.de
@*        Lars Kastner, lkastner@math.fu-berlin.de
@*        Oleksandr Motsak, U@D, where U={motsak}, D={mathematik.uni-kl.de}


OVERVIEW: using this library allows one can virtually add multigrading to Singular.
For more see http://code.google.com/p/convex-singular/wiki/Multigrading

NOTE: this is a proof of concept using Singular attributes. 

REFERENCE:
E. Miller, B. Sturmfels: 'Combinatorial Commutative Algebra', 
M. Kreuzer, L. Robbiano: 'Computational Commutative Algebra'

KEYWORDS:  multigradeding, multidegree, multiweights, multigraded-homogenous


REQUIRES: 4ti2 for Hilbert Bases needed for mDegBasis

PROCEDURES:
setBaseMultigrading(M,T) attach multiweights/torsion matrices to the basering
getVariableWeights([R])  get matrix of multidegrees of vars attached to a ring
getTorsion([R])          get torsion matrix attached to a ring

setModuleGrading(M,v)    attach multiweights of units to a module and return it
getModuleGrading(M)      get multiweights of module units (attached to M)

mDeg(A)                  compute the multidegree of A
mDegBasis(d)             compute all monomials of multidegree d
mDegPartition(p)         compute the multigraded-homogenous components of p

isTorsionFree()          test whether the current multigrading is torsion free
isTorsionElement(p)      test whether p has zero multidegree
isHomogenous(a)          test whether 'a' is multigraded-homogenous

equalMDeg(e1,e2[,V])     test whether e1==e2 in the current multigrading

mDegGroebner(M)          compute the multigraded GB/SB of M
mDegSyzygy(M)            compute the multigraded syzygies of M
mDegResolution(M,l[,m])  compute the multigraded resolution of M

defineHomogenous(p)      get a torsion matrix wrt which p becomes homogenous
pushForward(f)           find the finest grading on the image ring, homogenizing f

hermite(A)               compute the Hermite Normal Form of a matrix

hilbertSeries(M)         compute the multigraded Hilbert Series of M

           (parameters in square brackets [] are optional)
";

// finestMDeg(def r)
// newMap(map F, intmat Q, list #)



LIB "standard.lib"; // for groebner

/******************************************************/
proc setBaseMultigrading(intmat M, list #)
"USAGE:   setBaseMultigrading(M[, T])
PURPOSE: attach M abd T to the basering, where M is the matrix of the weights of variables, L is the torsion
RETURN: nothing
EXAMPLE: example setBaseMultigrading; shows an example
"
{
  string attrMgrad   = "mgrad";
  string attrTorsion = "torsion";
  string attrTorsionHNF = "torsion_hermite";


  attrib(basering, attrMgrad, M);
  
  if( size(#) > 0 and typeof(#[1]) == "intmat" )
  {
    def T = #[1];
  } else
  {
    intmat T[nrows(M)][1];
  }

  if( nrows(T) == nrows(M) )
  {
    attrib(basering, attrTorsion, T);  
    def H;
    attrib(basering, attrTorsionHNF, H);        
  }
  else
  {
    ERROR("Incompatible matrix sizes!");
  }
}
example
{
  "EXAMPLE:"; echo=2;
   
  ring R = 0, (x, y, z), dp;

  // Weights of variables
  intmat M[3][3] =
    1, 0, 0,
    0, 1, 0,
    0, 0, 1;

  // Torsion:
  intmat L[3][2] =
    1, 1,
    1, 3, 
    1, 5;
   
  // attaches M & L to R (==basering):
  setBaseMultigrading(M, L); // Grading: Z^3/L

  // Weights are accessible via "getVariableWeights()":
  getVariableWeights() == M;
  getVariableWeights(R) == M;
  getVariableWeights(basering) == M;

  // Torsion is accessible via "getTorsion()":
  getTorsion() == L;
  getTorsion(R) == L;
  getTorsion(basering) == L;

  // And its hermite NF via getTorsion("hermite"):
  intmat H = hermite(L);
  getTorsion("hermite") == H;
  getTorsion(R, "hermite") == H;
  getTorsion(basering, "hermite") == H;

  

  kill L, M;

  // ----------- isomorphic multigrading -------- //

  // Weights of variables
  intmat M[2][3] =
    1, -2, 1,
    1,  1, 0;

  // Torsion:
  intmat L[2][1] =
    0,
    2;
   
  // attaches M & L to R (==basering):
  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)

  // Weights are accessible via "getVariableWeights()":
  getVariableWeights() == M;

  // Torsion is accessible via "getTorsion()":
  getTorsion() == L;

  kill L, M;
  // ----------- extreme case ------------ //

  // Weights of variables
  intmat M[1][3] =
    1,  -1, 10;

  // Torsion:
  intmat L[1][1] =
    0;
   
  // attaches M & L to R (==basering):
  setBaseMultigrading(M); // Grading: Z^3

  // Weights are accessible via "getVariableWeights()":
  getVariableWeights() == M;

  // Torsion is accessible via "getTorsion()":
  getTorsion() == L;
}


/******************************************************/
proc getVariableWeights(list #)
"USAGE: getVariableWeights([R])
PURPOSE: get associated multigrading matrix for the basering [or R]
RETURN:  intmat, matrix of multidegrees of variables
EXAMPLE: example getVariableWeights; shows an example
"
{
  string attrMgrad = "mgrad";


  if( size(#) > 0 )
  {
    if(( typeof(#[1]) == "ring" ) || ( typeof(#[1]) == "qring" ))
    {
      def R = #[1];
    }
    else
    {
      ERROR("Optional argument must be a ring!");
    }
  }
  else
  {
    def R = basering;
  }

  def M = attrib(R, attrMgrad);
  if( typeof(M) == "intmat"){ return (M); }
  ERROR( "Sorry no multigrading matrix!" );  
}
example
{
  "EXAMPLE:"; echo=2;
   
  ring R = 0, (x, y, z), dp;

  // Weights of variables
  intmat M[3][3] =
    1, 0, 0,
    0, 1, 0,
    0, 0, 1;

  // Torsion:
  intmat L[3][2] =
    1, 1,
    1, 3, 
    1, 5;
   
  // attaches M & L to R (==basering):
  setBaseMultigrading(M, L); // Grading: Z^3/L

  // Weights are accessible via "getVariableWeights()":
  getVariableWeights() == M;

  kill L, M;

  // ----------- isomorphic multigrading -------- //

  // Weights of variables
  intmat M[2][3] =
    1, -2, 1,
    1,  1, 0;

  // Torsion:
  intmat L[2][1] =
    0,
    2;
   
  // attaches M & L to R (==basering):
  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)

  // Weights are accessible via "getVariableWeights()":
  getVariableWeights() == M;

  kill L, M;

  // ----------- extreme case ------------ //

  // Weights of variables
  intmat M[1][3] =
    1,  -1, 10;

  // Torsion:
  intmat L[1][1] =
    0;
   
  // attaches M & L to R (==basering):
  setBaseMultigrading(M); // Grading: Z^3

  // Weights are accessible via "getVariableWeights()":
  getVariableWeights() == M;
}

/******************************************************/
proc getTorsion(list #)
"USAGE: getTorsion([R[,opt]])
PURPOSE: get associated torsion matrix, i.e. generators (cols) of the Torsion group
RETURN: intmat, the torsion matrix, or its hermite normal form
        if an optional argument (\"hermite\") is given
"
{
  string attrTorsion    = "torsion";
  string attrTorsionHNF = "torsion_hermite";

  int i = 1;

  if( size(#) >= i )
  {
    if( ( typeof(#[i]) == "ring" ) or ( typeof(#[i]) == "qring" ) )
    {
      def R = #[i];
      i++;
    } 
  } 

  if( !defined(R) )
  {
    def R = basering;
  }

  if( size(#) >= i )
  {
    if( #[i] == "hermite" )
    {
      if( typeof(attrib(R, attrTorsionHNF)) != "intmat" ) 
      {
        def M = getTorsion(R);
        if( typeof(M) != "intmat")
        { 
          ERROR( "Sorry no torsion matrix!" );  
        }
        attrib(R, attrTorsionHNF, hermite(M)); // this might not work with R != basering...
      } 
      return (attrib(R, attrTorsionHNF));
    }
  }

  def M = attrib(R, attrTorsion);
  if( typeof(M) != "intmat")
  { 
    ERROR( "Sorry no torsion matrix!" );  
  }
  return (M);
}
example
{
  "EXAMPLE:"; echo=2;
   
  ring R = 0, (x, y, z), dp;

  // Weights of variables
  intmat M[3][3] =
    1, 0, 0,
    0, 1, 0,
    0, 0, 1;

  // Torsion:
  intmat L[3][2] =
    1, 1,
    1, 3, 
    1, 5;
   
  // attaches M & L to R (==basering):
  setBaseMultigrading(M, L); // Grading: Z^3/L

  // Torsion is accessible via "getTorsion()":
  getTorsion() == L;

  // its hermite NF:
  print(getTorsion("hermite"));

  kill L, M;

  // ----------- isomorphic multigrading -------- //

  // Weights of variables
  intmat M[2][3] =
    1, -2, 1,
    1,  1, 0;

  // Torsion:
  intmat L[2][1] =
    0,
    2;
   
  // attaches M & L to R (==basering):
  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)

  // Torsion is accessible via "getTorsion()":
  getTorsion() == L;

  // its hermite NF:
  print(getTorsion("hermite"));

  kill L, M;

  // ----------- extreme case ------------ //

  // Weights of variables
  intmat M[1][3] =
    1,  -1, 10;

  // Torsion:
  intmat L[1][1] =
    0;
   
  // attaches M & L to R (==basering):
  setBaseMultigrading(M); // Grading: Z^3

  // Torsion is accessible via "getTorsion()":
  getTorsion() == L;

  // its hermite NF:
  print(getTorsion("hermite"));
}

/******************************************************/
proc getModuleGrading(def m)
"USAGE: getModuleGrading(m), 'm' module/vector
RETURN: integer matrix of the multiweights of free module generators attached to 'm'
"
{
  string attrModuleGrading = "genWeights";

  //  print(m);  typeof(m);  attrib(m);

  def V = attrib(m, attrModuleGrading);

  if( typeof(V) != "intmat" )
  {
    if( (typeof(m) == "ideal") or (typeof(m) == "poly") )
    {
      intmat M = getVariableWeights();
      intmat VV[nrows(M)][1];
      return (VV);
    }
     
    ERROR("Sorry: vector or module need module-grading-matrix! See 'getModuleGrading'.");
  }

  if( nrows(V) != nrows(getVariableWeights()) )
  {
    ERROR("Sorry wrong height of V: " + string(nrows(V)));
  }

  if( ncols(V) < nrows(m) )
  {
    ERROR("Sorry wrong width of V: " + string(ncols(V)));
  }
   
  return (V);
}
example
{
  "EXAMPLE:"; echo=2;
   
   ring R = 0, (x,y), dp;
   intmat M[2][2]=
     1, 1,
     0, 2;
   intmat T[2][5]=
     1,  2,  3,  4, 0,
     0, 10, 20, 30, 1;
   
   setBaseMultigrading(M, T);
   
   ideal I = x, y, xy^5;
   isHomogenous(I);
 
   intmat V = mDeg(I); print(V);

   module S = syz(I); print(S);
   
   S = setModuleGrading(S, V);

   getModuleGrading(S) == V;
   
   vector v = setModuleGrading(S[1], V);
   getModuleGrading(v) == V;
   isHomogenous(v);
   print( mDeg(v) );   
   
   isHomogenous(S);
   print( mDeg(S) );
}

/******************************************************/
proc setModuleGrading(def m, intmat G)
"USAGE: setModuleGrading(m, G), m module/vector, G intmat
PURPOSE: attaches the multiweights of free module generators to 'm'
WARNING: The method does not verify that the multigrading makes the 
         module/vector homogenous. One can do that using isHomogenous(m).
"
{
  string attrModuleGrading = "genWeights";

  intmat R = getVariableWeights();

  if(nrows(G) != nrows(R)){ ERROR("Incompatible gradings.");}
  if(ncols(G) < nrows(m)){ ERROR("Multigrading does not fit to module.");}

  attrib(m, attrModuleGrading, G);
  return(m);
}
example
{
  "EXAMPLE:"; echo=2;
   
   ring R = 0, (x,y), dp;
   intmat M[2][2]=
     1, 1,
     0, 2;
   intmat T[2][5]=
     1,  2,  3,  4, 0,
     0, 10, 20, 30, 1;
   
   setBaseMultigrading(M, T);
   
   ideal I = x, y, xy^5;
   intmat V = mDeg(I);
   
   // V == M; modulo T
   print(V);

   module S = syz(I);
   
   S = setModuleGrading(S, V);
   getModuleGrading(S) == V;

   print(S);
   
   vector v = S[1]; v = setModuleGrading(v, V);
   getModuleGrading(v) == V;

   print( mDeg(v) );   

   isHomogenous(S);

   print( mDeg(S) );
}




/******************************************************/
proc isTorsionFree()
"USAGE: isTorsionFree()
PURPOSE: Determines whether the multigrading attached to the current ring is torsion-free.
RETURN: boolean, the result of the test
"
{

  intmat H = hermite(transpose(getTorsion("hermite"))); // TODO: ?cache it?

  int i, j;
  int r = nrows(H);
  int c = ncols(H);
  int d = 1;
  for( i = 1; (i <= c) && (i <= r); i++ )
  {
    for( j = i; (H[j, i] == 0)&&(j < r); j++ )
    {
    }

    if(H[j, i]!=0)
    { 
      d=d*H[j, i];
    }
  }

  if( (d*d)==1 )
  { 
    return(1==1);
  }
  return(0==1);
}
example
{
  "EXAMPLE:"; echo=2;

   ring R = 0,(x,y),dp;
   intmat M[2][2]=
     1,0,
     0,1;
   intmat T[2][5]=
     1, 2, 3, 4, 0,
     0,10,20,30, 1;
   
   setBaseMultigrading(M,T);
   
   // Is the resulting group torsion free?
   isTorsionFree();

   kill R, M, T;
   ///////////////////////////////////////////

   ring R=0,(x,y,z),dp;
   intmat A[3][3] = 
     1,0,0,
     0,1,0,
     0,0,1;
   intmat B[3][4]=
     3,3,3,3,
     2,1,3,0,
     1,2,0,3;
   setBaseMultigrading(A,B);
   // Is the resulting group torsion free?
   isTorsionFree();

   kill R, A, B;
}



/******************************************************/
proc hermite(intmat A)
"USAGE: hermite( A );
PROCEDURE: Computes the Hermite Normal Form of the matrix A by column operations.
RETURN: intmat, the Hermite Normal Form of A
"
{
  intmat g;
  int i,j,k, save, d, a1, a2, c, l;
  c=0;
  l=0;
  for(i=1; ((i+l)<=nrows(A))&&((i+l)<=ncols(A)); i++)
  {
    //i;
    if(A[i+l, i]==0)
    {
      for(j=i;j<=ncols(A); j++)
      {
        if(A[i+l, j]!=0)
        {
          for(k=1;k<=nrows(A);k++)
          {
            save=A[k, i];
            A[k, i]=A[k, j];
            A[k, j]=save;
          }
          break;
        }
        if(j==ncols(A)){ c=1; l=l+1; }
      }
    }

    if((i+l)>nrows(A)){ break; }
    
    if(A[i+l, i]<0)
    {
      for(k=1;k<=nrows(A);k++){ A[k, i]=A[k, i]*(-1); }
    }

    if(c==0)
    {
      for(j=i+1;j<=ncols(A);j++)
      {
        //A;
        if(A[i+l, j]<0)
        {
          for(k=1;k<=nrows(A);k++){ A[k, j]=(-1)*A[k, j];}
        }

        if(A[i+l, i]==1)
        {
          d=A[i+l, j];
          for(k=1;k<=nrows(A);k++)
          {
            A[k, j]=A[k, j]-d*A[k, i];
          }
        }
        else
        {
          while((A[i+l, i] * A[i+l, j])!=0)
          {
            if(A[i+l, i]> A[i+l, j])
            {

              for(k=1;k<=nrows(A);k++)
              { 
                save=A[k, i];
                A[k, i]=A[k, j];
                A[k, j]=save;
              }
            }
            a1=A[i+l, j]%A[i+l,i];
            a2=A[i+l, j]-a1;
            d=a2/A[i+l, i];
            for(k=1;k<=nrows(A);k++)
            { 
              A[k, j]=A[k, j]- d*A[k, i];
            }
          }
        }
      }
      for(j=1;j<i;j++)
      {
        a1=A[i+l, j]%A[i+l,i];
        d=(A[i+l, j]-a1)/A[i+l, i];
        for(k=1;k<=nrows(A);k++){ A[k, j]=A[k, j]-d*A[k, i];}
      }
    }
    c=0;
  }
  return( A);
}
example
{
  "EXAMPLE:"; echo=2;

   intmat M[2][5] = 
     1, 2, 3, 4, 0,
     0,10,20,30, 1;

   // Hermite Normal Form of M:
   print(hermite(M));

   intmat T[3][4] = 
     3,3,3,3,
     2,1,3,0,
     1,2,0,3;

   // Hermite Normal Form of T:
   print(hermite(T));

   intmat A[4][5] = 
     1,2,3,2,2,
     1,2,3,4,0,
     0,5,4,2,1,
     3,2,4,0,2;

   // Hermite Normal Form of A:
   print(hermite(A));
}


/******************************************************/
proc isTorsionElement(intvec mdeg)
"USAGE: isTorsionElement(intvec mdeg);
PROCEDURE: For a integer vector mdeg representing the multidegree of some polynomial
or vector this method computes if the multidegree is contained in the torsion
group, i.e. if it is zero in the multigrading.
"
{
  intmat H = getTorsion("hermite");
  int x, k, i;

  int r = nrows(H);
  int c = ncols(H);

  int rr = nrows(mdeg);

  for( i = 1; (i<=r) && (i<=c); i++)
  {
    if(H[i, i]!=0)
    {
      x = mdeg[i]%H[i, i];

      if(x!=0)
      { 
        return(1==0);
      }

      x = mdeg[i] / H[i,i];

      for( k=1; k <= rr; k++)
      { 
        mdeg[k] = mdeg[k] - x*H[k,i];
      }
    }
  }

  return( mdeg == 0 );

}
example
{
 "EXAMPLE:"; echo=2;

  ring r = 0,(x,y,z),dp;

  intmat g[2][3]=
    1,0,1,
    0,1,1;
  intmat t[2][1]=
    -2,
    1;

  setBaseMultigrading(g,t);

  poly a = x10yz;
  poly b = x8y2z;
  poly c = x4z2;
  poly d = y5;
  poly e = x2y2;
  poly f = z2;

  intvec v1 = mDeg(a) - mDeg(b);
  v1;
  isTorsionElement(v1);
  
  intvec v2 = mDeg(a) - mDeg(c);
  v2;
  isTorsionElement(v2);
  
  intvec v3 = mDeg(e) - mDeg(f);
  v3;
  isTorsionElement(v3);
  
  intvec v4 = mDeg(c) - mDeg(d);
  v4;
  isTorsionElement(v4);
}


/******************************************************/
proc defineHomogenous(poly f, list #)
"USAGE: defineHomogenous(f[, G]); polynomial f, integer matrix G
PURPOSE: Yields a matrix which has to be appended to the torsion matrix to make the 
polynomial f homogenous in the grading by grad.
"
{
  if( size(#) > 0 ) 
  {
    if( typeof(#[1]) == "intmat" ) 
    {
      intmat grad = #[1];
    }
  }

  if( !defined(grad) )
  {
    intmat grad = getVariableWeights();
  }

  intmat newtor[nrows(grad)][size(f)-1];
  int i,j;
  intvec l = grad*leadexp(f);
  intvec v;
  for(i=2; i <= size(f); i++)
  {
    v = grad * leadexp(f[i]) - l;
    for( j=1; j<=size(v); j++)
    {
      newtor[j,i-1] = v[j];
    }
  }
  return(newtor);
}
example
{
  "EXAMPLE:"; echo=2;

  ring r =0,(x,y,z),dp;
  intmat grad[2][3] = 
    1,0,1,
    0,1,1;

  setBaseMultigrading(grad);

  poly f = x2y3-z5+x-3zx;

  intmat M = defineHomogenous(f);
  M;
  defineHomogenous(f, grad) == M;
  
  isHomogenous(f);
  setBaseMultigrading(grad, M);
  isHomogenous(f);
}

/******************************************************/
proc pushForward(map f)
"USAGE: pushForward(f);
PURPOSE: Computes the finest grading of the image ring which makes the map f
a map of graded rings. The group map between the two grading groups is given
by transpose( (Id, 0) ). Pay attention that the group spanned by the columns of
the torsion matrix may not be a subgroup of the grading group. Still all columns
are needed to find the correct image of the preimage gradings.
"
{

  int k,i,j;
  f;

  intmat oldgrad=getVariableWeights(preimage(f));
  intmat oldtor=getTorsion(preimage(f));

  int n=nvars(preimage(f));
  int np=nvars(basering);
  int p=nrows(oldgrad);
  int pp=p+np;

  intmat newgrad[pp][np];

  for(i=1;i<=np;i++){ newgrad[p+i,i]=1;}

  //newgrad;



  list newtor;
  intmat toadd;
  int columns=0;

  intmat toadd1[pp][n];
  intvec v;
  poly im;

  for(i=1;i<=p;i++){
    for(j=1;j<=n;j++){ toadd1[i,j]=oldgrad[i,j];}
  }

  for(i=1;i<=n;i++){
    im=f[i];
    //im;
    toadd = defineHomogenous(im, newgrad);
    newtor=insert(newtor,toadd);
    columns=columns+ncols(toadd);

    v=leadexp(f[i]);
    for(j=p+1;j<=p+np;j++){ toadd1[j,i]=-v[j-p];}
  }

  newtor=insert(newtor,toadd1);
  columns=columns+ncols(toadd1);


  if(typeof(basering)=="qring"){
    //"Entering qring";
    ideal a=ideal(basering);
    for(i=1;i<=size(a);i++){
      toadd = defineHomogenous(a[i], newgrad);
      //toadd;
      columns=columns+ncols(toadd);
      newtor=insert(newtor,toadd);
    }
  }

  //newtor;
  intmat imofoldtor[pp][ncols(oldtor)];
  for(i=1; i<=nrows(oldtor);i++){
    for(j=1; j<=ncols(oldtor); j++){
      imofoldtor[i,j]=oldtor[i,j];
    }
  }

  columns=columns+ncols(oldtor);
  newtor=insert(newtor, imofoldtor);

  intmat torsion[pp][columns];
  columns=0;
  for(k=1;k<=size(newtor);k++){
    for(i=1;i<=pp;i++){
      for(j=1;j<=ncols(newtor[k]);j++){torsion[i,j+columns]=newtor[k][i,j];}
    }
    columns=columns+ncols(newtor[k]);
  }

  torsion=hermite(torsion);
  intmat result[pp][pp];
  for(i=1;i<=pp;i++){
    for(j=1;j<=pp;j++){result[i,j]=torsion[i,j];}
  }

  setBaseMultigrading(newgrad, result);

}
example
{
  "EXAMPLE:"; echo=2;

  // Setting degrees for preimage ring.;
  intmat grad[3][3] = 
    1,0,0,
    0,1,0,
    0,0,1;

  ring r = 0,(x,y,z),dp;
  keepring(r);
  setBaseMultigrading(grad);

  // grading on r:
  getVariableWeights();
  getTorsion();

  ring R = 0,(a,b),dp;
  ideal i=a2-b2+a6-b5+ab3,a7b+b15-ab6+a6b6;

  // The quotient ring by this ideal will become our image ring.;
  qring Q = std(i);
  map f = r,-a2b6+b5+a3b+a2+ab,-a2b7-3a2b5+b4+a,a6-b6-b3+a2;
  f;

  // Pushing forward f:
  pushForward(f);

  // due to pushForward we have got new grading on Q
  getVariableWeights();
  getTorsion();
  
  // TODO: Unfortunately this is not a very spectacular example.;

}


/******************************************************/
proc equalMDeg(intvec exp1, intvec exp2, list #)
"USAGE: equalMDeg(exp1, exp2[, V]) where exp1,exp2 are exponents of terms [, V is multidegree of module components];
PURPOSE: Tests if the exponent vectors of two monomials represent the same multidegree.
"
{
  if( size(exp1) != size(exp2) )
  {
    ERROR("Sorry: we cannot compare exponents comming from a polynomial and a vector yet!");
  }

  if( exp1 == exp2) 
  {
    return (1==1);
  }



  intmat M = getVariableWeights();

  if( nrows(exp1) > ncols(M) ) // vectors => last exponent is the module component!
  {
    if( (size(#) == 0) or (typeof(#[1])!="intmat") )
    {
      ERROR("Sorry: wrong or missing module-unit-weights-matrix!");
    }
    intmat V = #[1];

    // typeof(V); print(V);

    int N = ncols(M);
    int r = nrows(M);

    intvec d = intvec(exp1[1..N]) - intvec(exp2[1..N]);
    intvec dm = intvec(V[1..r, exp1[N+1]]) - intvec(V[1..r, exp2[N+1]]);

    intvec difference = M * d + dm;
  }
  else
  {
    intvec d = (exp1 - exp2);
    intvec difference = M * d;
  }

  if (isFreeRepresented()) // no torsion!?
  {
    return ( difference == 0);
  }
  return ( isTorsionElement( difference ) );
}
example
{
  "EXAMPLE:"; echo=2;

  ring r = 0,(x,y,z),dp;

  intmat g[2][3]=
    1,0,1,
    0,1,1;

  intmat t[2][1]=
    -2,
    1;

  setBaseMultigrading(g,t);

  poly a = x10yz;
  poly b = x8y2z;
  poly c = x4z2;
  poly d = y5;
  poly e = x2y2;
  poly f = z2;


  equalMDeg(leadexp(a), leadexp(b));
  equalMDeg(leadexp(a), leadexp(c));
  equalMDeg(leadexp(a), leadexp(d));
  equalMDeg(leadexp(a), leadexp(e));
  equalMDeg(leadexp(a), leadexp(f));

  equalMDeg(leadexp(b), leadexp(c));
  equalMDeg(leadexp(b), leadexp(d));
  equalMDeg(leadexp(b), leadexp(e));
  equalMDeg(leadexp(b), leadexp(f));

  equalMDeg(leadexp(c), leadexp(d));
  equalMDeg(leadexp(c), leadexp(e));
  equalMDeg(leadexp(c), leadexp(f));

  equalMDeg(leadexp(d), leadexp(e));
  equalMDeg(leadexp(d), leadexp(f));

  equalMDeg(leadexp(e), leadexp(f));

}



/******************************************************/
static proc isFreeRepresented()
"check whether the base muligrading is torsion-free (it is zero).
"
{
  intmat T = getTorsion();

  intmat Z[nrows(T)][ncols(T)];

  return (T == Z); // no torsion!
}


/******************************************************/
proc isHomogenous(def a, list #)
"USAGE: isHomogenous(a[, f]); a polynomial/vector/ideal/module
RETURN: boolean, TRUE if a is (multi)homogenous, and FALSE otherwise
"
{
  if( (typeof(a) == "poly") or (typeof(a) == "vector") )
  {
    return ( size(mDegPartition(a)) <= 1 )
  }

  if( typeof(a) == "vector" or typeof(a) == "module" )
  {
    intmat V = getModuleGrading(a);
  }

  if( (typeof(a) == "ideal") or (typeof(a) == "module") )
  {
    if(size(#) > 0)
    {
      if (#[1] == "checkGens")
      {
        def aa;
        for( int i = ncols(a); i > 0; i-- )
        {
          aa = a[i];

          if(defined(V))
          {
            aa = setModuleGrading(aa, V);
          }

          if(!isHomogenous(aa))
          {
            return(0==1);
          }
        }
        return(1==1);
      }
    }

    def g = groebner(a); // !!!!

    def b, aa; int j; 
    for( int i = ncols(a); i > 0; i-- )
    {
      aa = a[i];

      if(defined(V))
      {
        aa = setModuleGrading(aa, V);
      }

      b = mDegPartition(aa);
      for( j = ncols(b); j > 0; j-- )
      {
        if(NF(b[j],g) != 0)
        {
          return(0==1);
        }
      }
    }
    return(1==1);
  }  
}
example
{
  "EXAMPLE:"; echo=2;

  ring r = 0,(x,y,z),dp;

  //Grading and Torsion matrices:
  intmat M[3][3] = 
    1,0,0,
    0,1,0,
    0,0,1;

  intmat T[3][1] =
    1,2,3;

  setBaseMultigrading(M,T);

  attrib(r);

  poly f = x-yz;

  mDegPartition(f);
  print(mDeg(_));

  isHomogenous(f);   // f: is not homogenous

  poly g = 1-xy2z3;
  isHomogenous(g); // g: is homogenous
  mDegPartition(g);

  kill T;
  /////////////////////////////////////////////////////////
  // new Torsion matrix:
  intmat T[3][4] = 
    3,3,3,3,
    2,1,3,0,
    1,2,0,3;
  
  setBaseMultigrading(M,T);

  f;
  isHomogenous(f);
  mDegPartition(f);

  // --------------------- 
  g;
  isHomogenous(g);
  mDegPartition(g);

  kill r, T, M;

  ring R = 0, (x,y,z), dp;

  intmat A[2][3] = 
    0,0,1,
    3,2,1;
  intmat T[2][1] = 
    -1, 
     4;
  setBaseMultigrading(A, T);

  isHomogenous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3)); // 1
  isHomogenous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3), "checkGens");
  isHomogenous(ideal(x+y, x2 - y2)); // 0

  // Degree partition:
  mDegPartition(x2 - y3 -xy +z);
  mDegPartition(x3 -y2z + x2 -y3 + z + 1);

  
  module N = gen(1) + (x+y) * gen(2), z*gen(3);

  intmat V[2][3] = 0; // 1, 2, 3,  4, 5, 6; //  column-wise weights of components!!??
 
  vector v1, v2;
 
  v1 = setModuleGrading(N[1], V); v1;
  mDegPartition(v1);
  print( mDeg(_) );

  v2 = setModuleGrading(N[2], V); v2;
  mDegPartition(v2);
  print( mDeg(_) );

  N = setModuleGrading(N, V);
  isHomogenous(N);
  print( mDeg(N) );

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

  V = 
    1, 2, 3,  
    4, 5, 6;

  v1 = setModuleGrading(N[1], V); v1;
  mDegPartition(v1);
  print( mDeg(_) );

  v2 = setModuleGrading(N[2], V); v2;
  mDegPartition(v2);
  print( mDeg(_) );

  N = setModuleGrading(N, V);
  isHomogenous(N);
  print( mDeg(N) );

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

  V = 
    0, 0, 0,  
    4, 1, 0;

  N = gen(1) + x * gen(2), z*gen(3);

  N = setModuleGrading(N, V); N;
  isHomogenous(N);
  print( mDeg(N) );

  v1 = setModuleGrading(N[1], V); v1;
  mDegPartition(v1);
  print( mDeg(_) );

  N = setModuleGrading(N, V); N;
  isHomogenous(N);
  print( mDeg(N) );
}

/******************************************************/
proc mDeg(def A)
"USAGE: mDeg(A); any A
PURPOSE: compute multidegree
"
{
  if( defined(attrib(A, "grad")) > 0 )
  {
    return (attrib(A, "grad"));
  }

  intmat M = getVariableWeights();
  int N = nvars(basering);

  if( ncols(M) != N )
  {
    ERROR("Sorry wrong mgrad-size of M: " + string(ncols(M)));
  }

  int r = nrows(M);

  if( A == 0 )
  {
    intvec v; v[r] = 0;
    return (v);
  }

  if( (typeof(A) == "vector") or (typeof(A) == "module") )
  {
    intmat V = getModuleGrading(A);

    if( nrows(V) != r )
    {
      ERROR("Sorry wrong mgrad-size of V: " + string(nrows(V)));
    }
  }

  intvec m; m[r] = 0;

  if( typeof(A) == "poly" )
  {
    intvec v = leadexp(A); //  v;
    m = M * v;

    // We assume homogeneous input!
    return(m);

    A = A - lead(A);
    while( size(A) > 0 )
    { 
      v = leadexp(A); //  v;
      m = max( m, M * v, r ); // ????
      A = A - lead(A);
    }

    return(m);
  }


  if( typeof(A) == "vector" )
  {
    intvec v;
    v = leadexp(A); //  v;
    m = intvec(M * intvec(v[1..N])) + intvec(V[1..r, v[N+1]]);

    // We assume homogeneous input!
    return(m);

    A = A - lead(A);
    while( size(A) > 0 )
    { 
      v = leadexp(A); //  v;

      // intvec(M * intvec(v[1..N])) + intvec(V[1..r, v[N+1]]);

      m = max( m, intvec(M * intvec(v[1..N])) + intvec(V[1..r, v[N+1]]), r ); // ???

      A = A - lead(A);
    }

    return(m);
  }

  int i, j; intvec d;

  if( typeof(A) == "ideal" )
  {
    intmat G[ r ] [ ncols(A)];
    for( i = ncols(A); i > 0; i-- )
    {
      d = mDeg( A[i] );

      for( j = 1; j <= r; j++ ) // see ticket: 253
      {
        G[j, i] = d[j];
      }      
    }
    return(G);
  }

  if( typeof(A) == "module" )
  {
    intmat G[ r ] [ ncols(A)];
    vector v;

    for( i = ncols(A); i > 0; i-- )
    {
      v = setModuleGrading(A[i], V);

      // G[1..r, i] 
      d = mDeg(v);

      for( j = 1; j <= r; j++ ) // see ticket: 253
      {
        G[j, i] = d[j];
      }      

    }

    return(G);
  }
 
}
example
{
  "EXAMPLE:"; echo=2;

  ring r = 0,(x, y), dp;

  intmat A[2][2] = 1, 0, 0, 1;
  print(A);

  intmat Ta[2][1] = 0, 3;
  print(Ta);

  //   attrib(A, "torsion", Ta); // to think about

//  "poly:";
  setBaseMultigrading(A);


  mDeg( x*x, A );
  mDeg( y*y*y, A );

  setBaseMultigrading(A, Ta);
  
  mDeg( x*x*y );
  
  mDeg( y*y*y*x );
  
  mDeg( x*y + x + 1 );

  mDegPartition(x*y + x + 1);

  print ( mDeg(0) );
  poly zero = 0;
  print ( mDeg(zero) );

//  "ideal:";
  
  ideal I = y*x*x, x*y*y*y;
  print( mDeg(I) );

  print ( mDeg(ideal(0)) );
  print ( mDeg(ideal(0,0,0)) );

//  "vectors:";
  
  intmat B[2][2] = 0, 1, 1, 0;
  print(B);
  
  mDeg( setModuleGrading(y*y*y*gen(2), B ));
  mDeg( setModuleGrading(x*x*gen(1), B ));

  
  vector V = x*gen(1) + y*gen(2);
  V = setModuleGrading(V, B);
  mDeg( V );

  vector v1 = setModuleGrading([0, 0, 0], B);
  print( mDeg( v1 ) );
  
  vector v2 = setModuleGrading([0], B);
  print( mDeg( v2 ) );

//  "module:";
  
  module D = x*gen(1), y*gen(2);
  D;
  D = setModuleGrading(D, B);
  print( mDeg( D ) );
  

  module DD = [0, 0],[0, 0, 0];
  DD = setModuleGrading(DD, B);
  print( mDeg( DD ) );

  module DDD = [0, 0];
  DDD = setModuleGrading(DDD, B);
  print( mDeg( DDD ) );

};





/******************************************************/
proc mDegPartition(def p)
"USAGE: mDegPartition(def p), p polynomial/vector
RETURNS: an ideal/module consisting of multigraded-homogeneous parts of p
"
{
  if( typeof(p) == "poly" ) 
  {
    ideal I;      
    poly mp, t, tt;
  }
  else
  {
    if(  typeof(p) == "vector" )
    {
      module I;      
      vector mp, t, tt;
    }
    else
    {
      ERROR("Wrong ARGUMENT type!");
    }
  }

  if( typeof(p) == "vector" )
  {
    intmat V = getModuleGrading(p); 
  }
  else
  {
    intmat V;
  }

  if( size(p) > 1) 
  {
    intvec m;

    while( p != 0 )
    {
      m = leadexp(p);
      mp = lead(p); 
      p = p - lead(p);
      tt = p; t = 0;

      while( size(tt) > 0 )
      { 

        if( equalMDeg( leadexp(tt), m, V  ) ) // TODO: we make no caching of matrices (M,T,H,V), which remain the same!
        {
          mp = mp + lead(tt); // "mp", mp;
        }
        else
        {
          t = t + lead(tt);  //  "t", t;
        }

        tt = tt - lead(tt);
      }

      I[size(I)+1] = mp;

      p = t;
    }
  }
  else
  {
    I[1] = p; // single monom
  }

  if( typeof(I) == "module" )
  {
    I = setModuleGrading(I, V);
  }

  return (I);
}
example
{
  "EXAMPLE:"; echo=2;

  ring r = 0,(x,y,z),dp;

  intmat g[2][3]=
    1,0,1,
    0,1,1;
  intmat t[2][1]=
    -2,
    1;

  setBaseMultigrading(g,t);

  poly f = x10yz+x8y2z-x4z2+y5+x2y2-z2+x17z3-y6;

  mDegPartition(f);
  
  vector v = xy*gen(1)-x3y2*gen(2)+x4y*gen(3);
  intmat B[2][3]=1,-1,-2,0,0,1;
  v = setModuleGrading(v,B);
  getModuleGrading(v);
  
  mDegPartition(v, B);
}



/******************************************************/
static proc unitMatrix(int n)
{
  intmat A[n][n];
  
  for( int i = n; i > 0; i-- )
  {
    A[i,i] = 1;
  }

  return (A);
}



/******************************************************/
static proc finestMDeg(def r)
"
USAGE:
PURPOSE: finest multigrading
"
{
  def save = basering;
  setring (r);

  // in basering
      ideal I = ideal(basering);

  int n = 0; int i; poly p;
  for( i = ncols(I); i > 0; i-- )
  {
    p = I[i];
    if( size(p) > 1 )
    {
      n = n + (size(p) - 1);
    }
    else
    {
      I[i] = 0;
    }
  }

  int N = nvars(basering);
  intmat A = unitMatrix(N);



  if( n > 0) 
  {

    intmat L[N][n]; 
    //  list L;
    int j = n;

    for(  i = ncols(I); i > 0; i-- )
    {
      p = I[i];

      if( size(p) > 1 ) 
      {
        intvec m0 = leadexp(p);
        p = p - lead(p);

        while( size(p) > 0 )
        {
          L[ 1..N, j ] = leadexp(p) - m0;
          p = p - lead(p);
          j--;
        }
      }
    }

    print(L);
    setBaseMultigrading(A, L);      
  } 
  else
  {
    setBaseMultigrading(A);
  }

  //  ERROR("nope");

  //  ring T = integer, (x), (C, dp);

  setring(save);
  return (r);
}
example
{
  "EXAMPLE:"; echo=2;

  ring r = 0,(x, y), dp;
  qring q  = std(x^2 - y);

  finestMDeg(q);

}




/******************************************************/
static proc newMap(map F, intmat Q, list #)
"
USAGE: ?? no use now...
PURPOSE: ??
"
{
  attrib(F, "Q", Q);

  if( size(#) > 0 and typeof(#[1]) == "intmat" )
  { 
    attrib(F, "P", #[1]);
  }
  return (F);
}


/******************************************************/
static proc matrix2intmat( matrix M )
{
  execute( "intmat A[ "+ string(nrows(M)) + "]["+ string(ncols(M)) + "] = " + string(M) + ";" );
  return (A);
}


/******************************************************/
static proc leftKernelZ(intmat M)
"USAGE:  leftKernel(M);   M a matrix
RETURN:  module
PURPOSE: computes left kernel of matrix M (a module of all elements v such that vM=0)
EXAMPLE: example leftKernel; shows an example
"
{
  if( nameof(basering) != "basering" )
  {
    def save = basering;
  }

  ring r = integer, (x), dp;


  //  basering;
  module N = matrix((M)); // transpose
  //  print(N);

  def MM = modulo( N, std(0) ) ;
  //  print(MM);

  intmat R = (  matrix2intmat( MM ) ); // transpose

  if( defined(save) > 0 )
  {
    setring save;
  }

  kill r;
  return( R );
}
example
{
  "EXAMPLE:"; echo=2;

  ring r= 0,(x,y,z),dp;
  matrix M[3][1] = x,y,z;
  print(M);
  matrix L = leftKernel(M);
  print(L);
  // check:
  print(L*M);
};



/******************************************************/
// the following is taken from "sing4ti2.lib" as we need 'hilbert' from 4ti2

static proc hilbert4ti2intmat(intmat A, list #)
"USAGE:  hilbert4ti2(A[,i]);
@*       A=intmat
@*       i=int
ASSUME:  - A is a matrix with integer entries which describes the lattice
@*         as ker(A), if second argument is not present,
@*         as the left image Im(A) = {zA : z \in ZZ^k}, if second argument is a positive integer
@*       - number of variables of basering equals number of columns of A
@*         (for ker(A)) resp. of rows of A (for Im(A))
CREATE:  temporary files sing4ti2.mat, sing4ti2.lat, sing4ti2.mar
@*       in the current directory (I/O files for communication with 4ti2)
NOTE:    input rules for 4ti2 also apply to input to this procedure
@*       hence ker(A)={x|Ax=0} and Im(A)={xA}
RETURN:  toric ideal specified by Hilbert basis thereof
EXAMPLE: example graver4ti2; shows an example
"
{
//--------------------------------------------------------------------------
// Initialization and Sanity Checks
//--------------------------------------------------------------------------
   int i,j;
   int nr=nrows(A);
   int nc=ncols(A);
   string fileending="mat";
   if (size(#)!=0)
   {
//--- default behaviour: use ker(A) as lattice
//--- if #[1]!=0 use Im(A) as lattice
      if(typeof(#[1])!="int")
      {
         ERROR("optional parameter needs to be integer value");
      }
      if(#[1]!=0)
      {
         fileending="lat";
      }
   }
//--- we should also be checking whether all entries are indeed integers
//--- or whether there are fractions, but in this case the error message
//--- of 4ti2 is printed directly

//--------------------------------------------------------------------------
// preparing input file for 4ti2
//--------------------------------------------------------------------------
   link eing=":w sing4ti2."+fileending;
   string eingstring=string(nr)+" "+string(nc);
   write(eing,eingstring);
   for(i=1;i<=nr;i++)
   {
      kill eingstring;
      string eingstring;
      for(j=1;j<=nc;j++)
      {
        //          if(g(A[i,j])>0)||(char(basering)!=0)||(npars(basering)>0))
        //          {
        //             ERROR("Input to hilbert4ti2 needs to be a matrix with integer entries");
        //          }
        eingstring=eingstring+string(A[i,j])+" ";
      }
      write(eing, eingstring);
   }
   close(eing);

//----------------------------------------------------------------------
// calling 4ti2 and converting output
// Singular's string is too clumsy for this, hence we first prepare
// using standard unix commands
//----------------------------------------------------------------------
   j=system("sh","hilbert -q -n sing4ti2"); ////////// be quiet + no loggin!!!

   j=system("sh","awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.hil | sed s/[\\\ \\\t\\\v\\\f]/,/g | sed s/,+/,/g |sed s/,,/,/g|sed s/,,/,/g > sing4ti2.converted");
   if( defined(keepfiles) <= 0)
   {
      j=system("sh",("rm -f sing4ti2.hil sing4ti2."+fileending));
   }
//----------------------------------------------------------------------
// reading output of 4ti2
//----------------------------------------------------------------------
   link ausg=":r sing4ti2.converted";
//--- last entry ideal(0) is used to tie the list to the basering
//--- it will not be processed any further

   string s = read(ausg);
   string ergstr = "intvec erglist = " + s + "0;";
   execute(ergstr);
  
   //   print(erglist);
  
   int Rnc = erglist[1];
   int Rnr = erglist[2];
   
   intmat R[Rnr][Rnc];

   int k = 3;

   for(i=1;i<=Rnc;i++)
   {
     for(j=1;j<=Rnr;j++)
     {
       //       "i: ", i, ", j: ", j, ", v: ", erglist[k];
       R[j, i] = erglist[k];
       k = k + 1;
     }
   }

   return (R);
//--- get rid of leading entry 0;
//   toric=toric[2..ncols(toric)];
//   return(toric);
}
// A nice example here is the 3x3 Magic Squares
example
{
  "EXAMPLE:"; echo=2;

   ring r=0,(x1,x2,x3,x4,x5,x6,x7,x8,x9),dp;
   intmat M[7][9]=1,1,1,-1,-1,-1,0,0,0,1,1,1,0,0,0,-1,-1,-1,0,1,1,-1,0,0,-1,0,0,1,0,1,0,-1,0,0,-1,0,1,1,0,0,0,-1,0,0,-1,0,1,1,0,-1,0,0,0,-1,1,1,0,0,-1,0,-1,0,0;
   hilbert4ti2intmat(M);
   hermite(M);
}

/////////////////////////////////////////////////////////////////////////////
static proc getMonomByExponent(intvec exp)
{
  int n = nvars(basering);

  if( nrows(exp) < n )
  {
    n = nrows(exp);
  }

  poly m = 1; int e;

  for( int i = 1; i <= n; i++ )
  {
    e = exp[i];
    if( e < 0 )
    {
      ERROR("Negative exponent!!!");
    }

    m = m * (var(i)^e);
  }

  return (m);

}

/******************************************************/
proc mDegBasis(intvec d) 
"
USAGE: multidegree d
ASSUME: current ring is multigraded, monomial ordering is global
PURPOSE: compute all monomials of multidegree d
"
{
  def R = basering;  // setring R;

  intmat M = getVariableWeights(R);

  //  print(M);

  int nr = nrows(M);
  int nc = ncols(M);

  intmat A[nr][nc+1];
  A[1..nr, 1..nc] = M[1..nr, 1..nc];
  //typeof(A[1..nr, nc+1]);
  if( nr==1)
  {
    A[1..nr, nc+1]=-d[1];
  }
  else
  {
    A[1..nr, nc+1] = -d;
  }

  intmat T = getTorsion(R);

  if( isFreeRepresented() )
  {
    intmat B = hilbert4ti2intmat(A);

    //      matrix B = unitMatrix(nrows(T));
  }
  else
  {
    int n = ncols(T);

    nc = ncols(A);

    intmat AA[nr][nc + 2 * n];
    AA[1..nr, 1.. nc] = A[1..nr, 1.. nc]; 
    AA[1..nr, nc + (1.. n)] = T[1..nr, 1.. n]; 
    AA[1..nr, nc + n + (1.. n)] = -T[1..nr, 1.. n]; 


    //      print ( AA );

    intmat K = leftKernelZ(( AA ) ); // 

    //      print(K);

    intmat KK[nc][ncols(K)] = K[ 1.. nc, 1.. ncols(K) ];

    //      print(KK);
    //      "!";

    intmat B = hilbert4ti2intmat(transpose(KK), 1);  

    //      "!";      print(B);

  }


  //  print(A);



  int i;  
  int nnr = nrows(B);
  int nnc = ncols(B);
  ideal I, J;
  if(nnc==0){
    I=0;
    return(I);
  }
  I[nnc] = 0;
  J[nnc] = 0;

  for( i = 1; i <= nnc; i++ )
  {
    //      "i: ", i;    B[nnr, i];

    if( B[nnr, i] == 1)
    {
      // intvec(B[1..nnr-1, i]);
      I[i] = getMonomByExponent(intvec(B[1..nnr-1, i]));
    }
    else
    {
      if( B[nnr, i] == 0)
      {
        // intvec(B[1..nnr-1, i]);
        J[i] = getMonomByExponent(intvec(B[1..nnr-1, i]));
      }
    }
    //      I[i];
  }

  ideal Q = (ideal(basering));

  if ( size(Q) > 0 )
  {
    I = NF( I, lead(Q) );
    J = NF( J, lead(Q) ); // Global ordering!!!
  }

  I = simplify(I, 2); // d
  J = simplify(J, 2); // d

  attrib(I, "ZeroPart", J);

  return (I);

  //  setring ;
}
example
{
  "EXAMPLE:"; echo=2;

  ring R = 0, (x, y), dp;

  intmat g1[2][2]=1,0,0,1;
  intmat t[2][1]=2,0;
  intmat g2[2][2]=1,1,1,1;
  intvec v1=4,0;
  intvec v2=4,4;
  
  intmat g3[1][2]=1,1;
  setBaseMultigrading(g3);
  intvec v3=4:1;
  v3;
  mDegBasis(v3);
  
  setBaseMultigrading(g1,t);
  mDegBasis(v1);
  setBaseMultigrading(g2);
  mDegBasis(v2);
  
  intmat M[2][2] = 1, -1, -1, 1;
  intvec d = -2, 2;

  setBaseMultigrading(M);

  mDegBasis(d);
  attrib(_, "ZeroPart");

  kill R;
  ring R = 0, (x, y, z), dp;

  intmat M[2][3] = 1, -2, 1,     1, 1, 0;

  intmat T[2][1] = 0, 2;

  intvec d = 4, 1;

  setBaseMultigrading(M, T);

  mDegBasis(d);
  attrib(_, "ZeroPart");


  kill R;

  ring R = 0, (x, y, z), dp;
  qring Q = std(ideal( y^6+ x*y^3*z-x^2*z^2 ));


  intmat M[2][3] = 1, 1, 2,     2, 1, 1;
  //  intmat T[2][1] = 0, 2;

  setBaseMultigrading(M);

  intvec d = 6, 6;
  mDegBasis(d);
  attrib(_, "ZeroPart");



  kill R;
  ring R = 0, (x, y, z), dp;
  qring Q = std(ideal( x*z^3 - y *z^6, x*y*z  - x^4*y^2 ));


  intmat M[2][3] = 1, -2, 1,     1, 1, 0;
  intmat T[2][1] = 0, 2;

  intvec d = 4, 1;

  setBaseMultigrading(M, T);

  mDegBasis(d);
  attrib(_, "ZeroPart");
}


proc mDegSyzygy(def I)
"USAGE: mDegSyzygy(I); I is a poly/vector/ideal/module
PURPOSE: computes the multigraded syzygy of I
RETURNS: module, the syzygy of I
NOTE: generators of I must be multigraded homogeneous 
"
{
  if( isHomogenous(I, "checkGens") == 0) 
  { 
    ERROR ("Sorry: inhomogenous input!"); 
  } 
  module S = syz(I);
  S = setModuleGrading(S, mDeg(I));
  return (S);
}
example
{
  "EXAMPLE:"; echo=2;
  

  ring r = 0,(x,y,z,w),dp;

  intmat M[2][4]=
    1,1,1,1,
    0,1,3,4;

  setBaseMultigrading(M);

  
  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
  
 
  intmat v[2][nrows(M)]=
    1,
    0;
  
  M = setModuleGrading(M, v);

  isHomogenous(M);
  "Multidegrees: "; print(mDeg(M));

  // Let's compute Syzygy!
  def S = mDegSyzygy(M); S;
  "Module Units Multigrading: "; print( getModuleGrading(S) );
  "Multidegrees: "; print(mDeg(S));

  isHomogenous(S);
}


proc mDegGroebner(def I)
"USAGE: mDegGroebner(I); I is a poly/vector/ideal/module
PURPOSE: computes the multigraded standard/groebner basis of I
NOTE: I must be multigraded homogeneous 
RETURNS: ideal/module, the computed basis
"
{
  if( isHomogenous(I) == 0) 
  { 
    ERROR ("Sorry: inhomogenous input!"); 
  } 

  def S = groebner(I);
  
  if( typeof(I) == "module" or typeof(I) == "vector" )
  {
    S = setModuleGrading(S, getModuleGrading(I));      
  }

  return(S);
}
example
{
  "EXAMPLE:"; echo=2;

  ring r = 0,(x,y,z,w),dp;

  intmat M[2][4]=
    1,1,1,1,
    0,1,3,4;

  setBaseMultigrading(M);

  
  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
  
 
  intmat v[2][nrows(M)]=
    1,
    0;
  
  M = setModuleGrading(M, v);


  /////////////////////////////////////////////////////////////////////////////
  // GB:
  M = mDegGroebner(M); M;
  "Module Units Multigrading: "; print( getModuleGrading(M) );
  "Multidegrees: "; print(mDeg(M));

  isHomogenous(M);

  /////////////////////////////////////////////////////////////////////////////
  // Let's compute Syzygy!
  def S = mDegSyzygy(M); S;
  "Module Units Multigrading: "; print( getModuleGrading(S) );
  "Multidegrees: "; print(mDeg(S));

  isHomogenous(S);

  /////////////////////////////////////////////////////////////////////////////
  // GB:
  S = mDegGroebner(S); S;
  "Module Units Multigrading: "; print( getModuleGrading(S) );
  "Multidegrees: "; print(mDeg(S));

  isHomogenous(S);
}


/******************************************************/
proc mDegResolution(def I, int ll, list #)
"USAGE: mDegResolution(I,l,[f]); I is poly/vector/ideal/module; l,f are integers
PURPOSE: computes the multigraded resolution of I of the length l, 
or the whole resolution if l is zero. Returns minimal resolution if an optional 
argument 1 is supplied
NOTE: input must have multigraded-homogeneous generators.
The returned list is trunkated at the first zero entry.
RETURNS: list, the computed resolution
"
{
  if( isHomogenous(I, "checkGens") == 0) 
  { 
    ERROR ("Sorry: inhomogenous input!"); 
  } 

  def R = res(I, ll, #); list L = R; int l = size(L);

  if( (typeof(I) == "module") or (typeof(I) == "vector") )
  {
    L[1] = setModuleGrading(L[1], getModuleGrading(I));      
  }

  int i; 
  for( i = 2; i <= l; i++ )
  {
    if( size(L[i]) > 0 )
    {
      L[i] = setModuleGrading( L[i], mDeg(L[i-1]) );
    } else
    {
      return (L[1..(i-1)]);
    }
  }
 
  return (L);

  
}
example
{
  "EXAMPLE:"; echo=2;
  
  ring r = 0,(x,y,z,w),dp;

  intmat M[2][4]=
    1,1,1,1,
    0,1,3,4;

  setBaseMultigrading(M);

  
  module m= ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
  
  isHomogenous(ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3), "checkGens");
  
  ideal A = xw-yz, x2z-y3, xz2-y2w, yw2-z3;

  int j;
  
  for(j=1; j<=ncols(A); j++)
  {
    mDegPartition(A[j]);
  }
  
  intmat v[2][1]=
    1,
    0;
  
  m = setModuleGrading(m, v);

  // Let's compute Syzygy!
  def S = mDegSyzygy(m); S;
  "Module Units Multigrading: "; print( getModuleGrading(S) );
  "Multidegrees: "; print(mDeg(S));

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

  S = mDegGroebner(S); S;
  "Module Units Multigrading: "; print( getModuleGrading(S) );
  "Multidegrees: "; print(mDeg(S));

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

  def L = mDegResolution(m, 0, 1);

  for( j =1; j<=size(L); j++)
  {
    "----------------------------------- ", j, " -----------------------------";
    L[j];
    "Module Multigrading: "; print( getModuleGrading(L[j]) );
    "Multigrading: "; print(mDeg(L[j]));
  }

  /////////////////////////////////////////////////////////////////////////////
  
  L = mDegResolution(maxideal(1), 0, 1);

  for( j =1; j<=size(L); j++)
  {
    "----------------------------------- ", j, " -----------------------------";
    L[j];
    "Module Multigrading: "; print( getModuleGrading(L[j]) );
    "Multigrading: "; print(mDeg(L[j]));
  }
  
  kill v;
  

  def h = hilbertSeries(m);
  setring h;

  numerator1;
  factorize(numerator1);
  
  denominator1;
  factorize(denominator1);

  numerator2;
  factorize(numerator2);

  denominator2;
  factorize(denominator2);
}

/******************************************************/
proc hilbertSeries(def I)
"USAGE: hilbertSeries(I); I is poly/vector/ideal/module
PURPOSE: computes the multigraded Hilbert Series of M 
NOTE: input must have multigraded-homogeneous generators. 
Multigrading should be positive.
RETURNS: a ring in variables t_(i), s_(i), with polynomials 
numerator1 and denominator1 and muturally prime numerator2 
and denominator2, quotients of which give the series.
"
{
   
  if( !isFreeRepresented() )
  {
    ERROR("SORRY: ONLY TORSION-FREE CASE (POSITIVE GRADING)");
  }
   
  int i, j, k, v;

  intmat M = getVariableWeights();
  
  int cc = ncols(M);
  int n = nrows(M);

  if( n == 0 )
  {
    ERROR("Error: wrong Variable Weights?");
  }

  list RES = mDegResolution(I,0,1);

  int l = size(RES);
  
  list L; L[l + 1] = 0;

  if(typeof(I) == "ideal")
  {
    intmat zeros[n][1];
    L[1] = zeros;
  } 
  else
  {
    L[1] = getModuleGrading(RES[1]);
  }

  for( j = 1; j <= l; j++)
  {
    L[j + 1] = mDeg(RES[j]);
  }
  
  l++;

  ring R = 0,(t_(1..n),s_(1..n)),dp;
  
  ideal units; 
  for( i=n; i>=1; i--)
  {
    units[i] = (var(i) * var(n + i) - 1);
  }
  
  qring Q = std(units);
  
  // TODO: should not it be a quotient ring depending on Torsion???
  // I am not sure about what to do in the torsion case, but since
  // we want to evaluate the polynomial at certain points to get
  // a dimension we need uniqueness for this. I think we would lose
  // this uniqueness if switching to this torsion ring.

  poly monom, summand, numerator;
  poly denominator = 1;
  
  for( i = 1; i <= cc; i++)
  {
    monom = 1;

    for( k = 1; k <= n; k++)
    {
      v = M[k,i];

      if(v >= 0)
      {
        monom = monom * (var(k)^(v));
      } 
      else
      {
        monom = monom * (var(n+k)^(-v));
      }
    }
    
    if( monom == 1)
    {
      ERROR("Multigrading not positive.");
    }

    denominator = denominator * (1 - monom);
  }
  
  for( j = 1; j<= l; j++) 
  {
    summand = 0;
    M = L[j];

    for( i = 1; i <= ncols(M); i++)
    {
      monom = 1;
      for( k = 1; k <= n; k++)
      {
        v = M[k,i];
        if( v > 0 )
        {
          monom = monom * (var(k)^v);
        } 
        else
        {
          monom = monom * (var(n+k)^(-v));
        }
      }
      summand = summand + monom;
    }
    numerator = numerator - (-1)^j * summand;
  }
  
  if( denominator == 0 )
  {
    ERROR("Multigrading not positive.");
  } 
  
  poly denominator1 = denominator;
  poly numerator1 = numerator;

  export denominator1;
  export numerator1;

  if( numerator != 0 )
  {
    poly d = gcd(denominator, numerator);

    poly denominator2 = denominator/d;
    poly numerator2 = numerator/d;

    if( gcd(denominator2, numerator2) != 1 )
    {
      ERROR("Sorry: gcd should be 1 (after dividing out gcd)! Something went wrong!");
    }
  }
  else
  {
    poly denominator2 = denominator;
    poly numerator2 = numerator;
  }


  export denominator2;
  export numerator2;

  " ------------ ";
  "This proc returns a ring with polynomials called 'numerator1/2' and 'denominator1/2'!";
  "They represent the first and the second Hilbert Series.";
  "The s_(i)-variables are defined to be the inverse of the t_(i)-variables.";
  " ------------ ";
  
  return(Q);
}
example
{
  "EXAMPLE:"; echo=2;
  
  ring r = 0,(x,y,z,w),dp;
  intmat g[2][4]=
    1,1,1,1,
    0,1,3,4;
  setBaseMultigrading(g);
  
  module M = ideal(xw-yz, x2z-y3, xz2-y2w, yw2-z3);
  intmat V[2][1]=
    1,
    0;

  M = setModuleGrading(M, V);

  def h = hilbertSeries(M); setring h;

  factorize(numerator2);
  factorize(denominator2);
  
  kill g, h; setring r;

  intmat g[2][4]=
    1,2,3,4,
    0,0,5,8;
  
  setBaseMultigrading(g);
  
  ideal I = x^2, y, z^3;
  I = std(I);
  def L = mDegResolution(I, 0, 1);

  for( int j = 1; j<=size(L); j++)
  {
    "----------------------------------- ", j, " -----------------------------";
    L[j];
    "Module Multigrading: "; print( getModuleGrading(L[j]) );
    "Multigrading: "; print(mDeg(L[j]));
  }

  mDeg(I);
  def h = hilbertSeries(I); setring h;
  
  factorize(numerator2);
  factorize(denominator2);

  kill r, h, g, V;
  ////////////////////////////////////////////////
  ////////////////////////////////////////////////
  
  ring R = 0,(x,y,z),dp;
  intmat W[2][3] = 
     1,1, 1,
     0,0,-1;
  setBaseMultigrading(W);
  ideal I = x3y,yz2,y2z,z4;
  
  def h = hilbertSeries(I); setring h;
  
  factorize(numerator2);
  factorize(denominator2);

  kill R, W, h;
  ////////////////////////////////////////////////
  ////////////////////////////////////////////////

  ring R = 0,(x,y,z,a,b,c),dp;
  intmat W[2][6] = 
     1,1, 1,1,1,1,
     0,0,-1,0,0,0;
  setBaseMultigrading(W);
  ideal I = x3y,yz2,y2z,z4;
  
  def h = hilbertSeries(I); setring h;
  
  factorize(numerator2);
  factorize(denominator2);
  
  kill R, W, h;

  ////////////////////////////////////////////////
  ////////////////////////////////////////////////
  ////////////////////////////////////////////////
  // This is example 5.3.9. from Robbianos book.
  
  ring R = 0,(x,y,z,w),dp;
  intmat W[1][4] = 
     1,1, 1,1;
  setBaseMultigrading(W);
  ideal I = z3,y3zw2,x2y4w2xyz2;

  hilb(std(I));
  
  def h = hilbertSeries(I); setring h;
  
  numerator1;
  denominator1;

  factorize(numerator2);
  factorize(denominator2);
  

  kill h;
  ////////////////////////////////////////////////
  setring R;

  ideal I2 = x2,y2,z2; I2;

  hilb(std(I2));
  
  def h = hilbertSeries(I2); setring h;

  numerator1;
  denominator1;


  kill h;
  ////////////////////////////////////////////////
  setring R;
  
  W = 2,2,2,2;
  
  setBaseMultigrading(W);

  getVariableWeights();

  intvec w = 2,2,2,2;

  hilb(std(I2), 1, w);

  kill w;
  

  def h = hilbertSeries(I2); setring h;

  
  numerator1; denominator1;
  kill h;

  
  kill R, W;

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

  ring R = 0,(x),dp;
  intmat W[1][1] =
     1;
  setBaseMultigrading(W);

  ideal I;

  I = 1; I;

  hilb(std(I));
  
  def h = hilbertSeries(I); setring h;

  numerator1; denominator1;

  kill h;
  ////////////////////////////////////////////////
  setring R;

  I = x; I;

  hilb(std(I));

  def h = hilbertSeries(I); setring h;

  numerator1; denominator1;
  
  kill h;  
  ////////////////////////////////////////////////
  setring R;

  I = x^5; I;

  hilb(std(I));
  hilb(std(I), 1);

  def h = hilbertSeries(I); setring h;

  numerator1; denominator1;
  
  
  kill h;  
  ////////////////////////////////////////////////
  setring R;

  I = x^10; I;

  hilb(std(I));

  def h = hilbertSeries(I); setring h;

  numerator1; denominator1;

  kill h;
  ////////////////////////////////////////////////
  setring R;

  module M = 1;

  M = setModuleGrading(M, W);

  
  hilb(std(M));

  def h = hilbertSeries(M); setring h;

  numerator1; denominator1;

  kill h;
  ////////////////////////////////////////////////
  setring R;

  kill M; module M = x^5*gen(1);

//  intmat V[1][3] = 0; // TODO: this would lead to a wrong result!!!?

  intmat V[1][1] = 0; // all gen(i) of degree 0!

  M = setModuleGrading(M, V);

  hilb(std(M));

  def h = hilbertSeries(M); setring h;

  numerator1; denominator1;

  kill h;
  ////////////////////////////////////////////////
  setring R; 

  module N = x^5*gen(3);

  kill V;
  
  intmat V[1][3] = 0; // all gen(i) of degree 0!

  N = setModuleGrading(N, V);
      
  hilb(std(N));

  def h = hilbertSeries(N); setring h;

  numerator1; denominator1;

  kill h;
  ////////////////////////////////////////////////
  setring R; 

  
  module S = M + N;
  
  S = setModuleGrading(S, V);

  hilb(std(S));

  def h = hilbertSeries(S); setring h;

  numerator1; denominator1;

  kill h;

  kill V;
  kill R, W;

}



///////////////////////////////////////////////////////////////////////////////
// testing for consistency of the library:
proc testMultigradingLib ()
{
  example setBaseMultigrading;
  example setModuleGrading;

  example getVariableWeights;
  example getTorsion;
  example getModuleGrading;


  example mDeg;
  example mDegPartition;


  example hermite;
  example isHomogenous;
  example isTorsionFree;
  example pushForward;
  example defineHomogenous;

  example equalMDeg;
  example isTorsionElement;

  example mDegResolution;
  example hilbertSeries;


// example mDegBasis; // needs 4ti2!
}