Index: Singular/ChangeLog
===================================================================
--- Singular/ChangeLog	(revision 807759680bf2373d9851bbe0a7f175e190f4c18a)
+++ Singular/ChangeLog	(revision 0879467b3486695745e782f3a4ad5ef36e73127d)
@@ -1,5 +1,5 @@
 - new command: rank
 - function: datetime()
-- new libs: resjung.lib, monomial.lib, surfacesignature.lib
+- new libs: resjung.lib, monomial.lib, surfacesignature.lib, multigrading.lib
 
 Fr 20. Aug 13:36:49 CEST 2010
Index: Singular/LIB/COPYING
===================================================================
--- Singular/LIB/COPYING	(revision 807759680bf2373d9851bbe0a7f175e190f4c18a)
+++ Singular/LIB/COPYING	(revision 0879467b3486695745e782f3a4ad5ef36e73127d)
@@ -114,6 +114,9 @@
                 Ph.Gimenez                      pgimenez@agt.uva.es
                 G.-M.Greuel                     greuel@mathematik.uni-kl.de
+multigrading.lib Rene Birkner                   rbirkner@math.fu-berlin.de
+                Lars Kastner                    lkastner@math.fu-berlin.de
+                Oleksandr Motsak                motsak at mathematik dot uni minus kl dot de
 ncalg.lib       Viktor Levandovskyy             levandov@mathematik.uni-kl.de
-                Oleksandr Motsak                motsak@mathematik.uni-kl.de
+                Oleksandr Motsak                motsak at mathematik dot uni minus kl dot de
 ncdecomp.lib    Viktor Levandovskyy             levandov@mathematik.uni-kl.de
 nchomolog.lib   Viktor Levandovskyy             levandov@math.rwth-aachen.de
@@ -122,5 +125,5 @@
                 Lobillo, F.J.                   jlobillo@ugr.es
                 Rabelo, C.                      crabelo@ugr.es
-                Oleksandr Motsak                motsak@mathematik.uni-kl.de
+                Oleksandr Motsak                motsak at mathematik dot uni minus kl dot de
 noether.lib     A. Hashemi                      Amir.Hashemi@lip6.fr
 normal.lib      Gert-Martin Greuel              greuel@mathematik.uni-kl.de
@@ -131,5 +134,5 @@
                 Dietmar Hillebrand
 paramet.lib     Thomas Keilen                   keilen@mathematik.uni-kl.de
-perron.lib      Oleksandr Motsak                motsak@mathematik.uni-kl.de
+perron.lib      Oleksandr Motsak                motsak at mathematik dot uni minus kl dot de
 phindex.lib 
 pointid.lib     Stefan Steidel                  steidel@mathematik.uni-kl.de
Index: Singular/LIB/multigrading.lib
===================================================================
--- Singular/LIB/multigrading.lib	(revision 0879467b3486695745e782f3a4ad5ef36e73127d)
+++ Singular/LIB/multigrading.lib	(revision 0879467b3486695745e782f3a4ad5ef36e73127d)
@@ -0,0 +1,2843 @@
+///////////////////////////////////////////////////////////////////
+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!
+}
Index: Singular/singular-libs
===================================================================
--- Singular/singular-libs	(revision 807759680bf2373d9851bbe0a7f175e190f4c18a)
+++ Singular/singular-libs	(revision 0879467b3486695745e782f3a4ad5ef36e73127d)
@@ -33,5 +33,5 @@
 # the will be included in share.tar.gz, but not in all.lib
 
-SLIBS1 = findiff.lib resjung.lib surfacesignature.lib
+SLIBS1 = findiff.lib resjung.lib surfacesignature.lib multigrading.lib
 
 PLIBS = bfun.lib  central.lib  dmod.lib  dmodapp.lib dmodvar.lib \
Index: doc/NEWS.texi
===================================================================
--- doc/NEWS.texi	(revision 807759680bf2373d9851bbe0a7f175e190f4c18a)
+++ doc/NEWS.texi	(revision 0879467b3486695745e782f3a4ad5ef36e73127d)
@@ -63,4 +63,6 @@
 @c @item
 @c new library: findiff.lib (see @nref{findiff_lib}: finite difference schemes)
+@c @item
+@c new library: multigrading.lib (see @nref{multigrading_lib}: multigradings and related computations)
 @end itemize
 
