///////////////////////////////////////////////////////////////////
version="$Id$";
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: This library allows one to virtually add multigradings to Singular.
For more see http://code.google.com/p/convex-singular/wiki/Multigrading
For theoretical references see:
E. Miller, B. Sturmfels: 'Combinatorial Commutative Algebra' and
M. Kreuzer, L. Robbiano: 'Computational Commutative Algebra'.

NOTE: 'mDegBasis' relies on 4ti2 for computing Hilbert Bases.

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)

KEYWORDS:  multigradeding, multidegree, multiweights, multigraded-homogenous
";

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

LIB "standard.lib"; // for groebner

/******************************************************/
proc setBaseMultigrading(intmat M, list #)
"USAGE: setBaseMultigrading(M[, T]); M, T are integer matrices
PURPOSE: attaches weights of variables and torsion to the basering.
NOTE: M encodes the weights of variables column-wise.
The torsion is given by the lattice spanned by the columns of the integer
matrix T in Z^nrows(M) over Z.
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();

  // Test all possible usages:
  (getVariableWeights() == M) && (getVariableWeights(R) == M) && (getVariableWeights(basering) == M);

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

  // Test all possible usages:
  (getTorsion() == L) && (getTorsion(R) == L) && (getTorsion(basering) == L);

  // And its hermite NF via getTorsion("hermite"):
  getTorsion("hermite");

  // Test all possible usages:
  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
EXAMPLE: example getTorsion; shows an example
"
{
  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!" );
        }
        M = hermite(M);
        attrib(R, attrTorsionHNF, 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'
EXAMPLE: example getModuleGrading; shows an example
"
{
  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 whether the multigrading makes the
         module/vector homogenous. One can do that using isHomogenous(m).
EXAMPLE: example setModuleGrading; shows an example
"
{
  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
EXAMPLE: example isTorsionFree; shows an example
"
{

  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 );
PURPOSE: Computes the (lower triangular) Hermite Normal Form
           of the matrix A by column operations.
RETURN: intmat, the Hermite Normal Form of A
EXAMPLE: example hermite; shows an example
"
{
  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);
PURPOSE: 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.
EXAMPLE: example isTorsionElement; shows an example
"
{
  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.
EXAMPLE: example defineHomogenous; shows an example
"
{
  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.
EXAMPLE: example pushForward; shows an example
"
{

  int k,i,j;
//  f;

//  listvar();
  def pre = preimage(f);

//  "pre: ";  pre;

  intmat oldgrad=getVariableWeights(pre);
  intmat oldtor=getTorsion(pre);

  int n=nvars(pre);
  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;

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



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

  setBaseMultigrading(grad);

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

  // only for the purpose of this example
  if( voice > 1 ){ /*keepring(r);*/ export(r); }

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

  listvar();

  map f = r,-a2b6+b5+a3b+a2+ab,-a2b7-3a2b5+b4+a,a6-b6-b3+a2; f;


  // TODO: Unfortunately this is not a very spectacular example...:
  // Pushing forward f:
  pushForward(f);

  // due to pushForward we have got new grading on Q
  getVariableWeights();
  getTorsion();


  // only for the purpose of this example
  if( voice > 1 ){ kill r; }

}


/******************************************************/
proc equalMDeg(intvec exp1, intvec exp2, list #)
"USAGE: equalMDeg(exp1, exp2[, V]); intvec exp1, exp2, intmat V
PURPOSE: Tests if the exponent vectors of two monomials (given by exp1 and exp2)
represent the same multidegree.
NOTE: the integer matrix V encodes multidegrees of module components,
if module component is present in exp1 and exp2
EXAMPLE: example equalMDeg; shows an example
"
{
  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 V!");
    }
    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
EXAMPLE: example isHomogenous; shows an example
"
{
  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); print(N);
  isHomogenous(N);
  print( mDeg(N) );
  v1 = setModuleGrading(N[1], V); print(v1);
  mDegPartition(v1);
  print( mDeg(_) );
  N = setModuleGrading(N, V); print(N);
  isHomogenous(N);
  print( mDeg(N) );
}

/******************************************************/
proc mDeg(def A)
"USAGE: mDeg(A); polynomial/vector/ideal/module A
PURPOSE: compute multidegree
EXAMPLE: example mDeg; shows an example
"
{
  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
EXAMPLE: example mDegPartition; shows an example
"
{
  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 )
      {
        // TODO: we do not cache matrices (M,T,H,V), which remain the same :(
        // TODO: we need some low-level procedure with all these arguments...!
        if( equalMDeg( leadexp(tt), m, V  ) )
        {
          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: finestMDeg(r); ring r
RETURN: ring, r endowed with the finest multigrading
TODO: not yet...
"
{
  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: newMap(F, Q[, P]); map F, intmat Q[, intmat P]
PURPOSE: endowe the map F with the integer matrices P [and Q]
"
{
  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, and
@*         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
EXAMPLE: example mDegBasis; shows an example
"
{
  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 module of I
RETURNS: module, the syzygy module of I
NOTE: generators of I must be multigraded homogeneous
EXAMPLE: example mDegSyzygy; shows an example
"
{
  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 MM[2][4]=
    1,1,1,1,
    0,1,3,4;
  setBaseMultigrading(MM);
  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 syzygies!
  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
EXAMPLE: example mDegGroebner; shows an example
"
{
  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 MM[2][4]=
    1,1,1,1,
    0,1,3,4;

  setBaseMultigrading(MM);


  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 truncated beginning with the first zero differential.
RETURNS: list, the computed resolution
EXAMPLE: example mDegResolution; shows an example
"
{
  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.
EXAMPLE: example hilbertSeries; shows an example
"
{

  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 ()
{
  echo = 2; printlevel = 3;
  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 hilbertSeries;


// example mDegBasis; // needs 4ti2!

  "The End!";

}