Changeset b6ae8c in git

Timestamp:

Mar 15, 2011, 7:17:52 PM (12 years ago)

Author:

Oleksandr Motsak <motsak@…>

Branches:

(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')

Children:

199cd68d23f96331d3327f7a0b11a46accefa1db

Parents:

4f7d761d43192f9c176ef0f44b0f61e891545eb8

Message:

ADD: newest version of multigrading.lib

From: Oleksandr Motsak <motsak@mathematik.uni-kl.de>

git-svn-id: file:///usr/local/Singular/svn/trunk@13984 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• Singular/LIB/multigrading.lib (modified) (147 diffs)

Singular/LIB/multigrading.lib

@@ +1 @@
+// -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// Vi-modeline: vim: filetype=c:syntax:shiftwidth=2:tabstop=8:textwidth=0:expandtab
 ///////////////////////////////////////////////////////////////////
 version="$Id$";
@@ -5 +7 @@
 LIBRARY:  multigrading.lib          Multigraded Rings
 
-AUTHORS:  Rene Birkner, rbirkner@math.fu-berlin.de
+AUTHORS:  Benjamin Bechtold, benjamin.bechtold@googlemail.com
+@*        Rene Birkner, rbirkner@math.fu-berlin.de
 @*        Lars Kastner, lkastner@math.fu-berlin.de
+@*        Simon Keicher, keicher@mail.mathematik.uni-tuebingen.de
 @*        Oleksandr Motsak, U@D, where U={motsak}, D={mathematik.uni-kl.de}
 
@@ -12 +16 @@
 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:
+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.
+All groups are finitely generated Abelian
 
 PROCEDURES:
-setBaseMultigrading(M,T); attach multiweights/torsion matrices to the basering
+setBaseMultigrading(M,L); attach multiweights/grading group 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
+
+getGradingGroup([R]);     get grading group attached to a ring
+getLattice([R[,choice]]); get grading group' lattice attached to a ring (or its NF)
+
+createGroup(S,L);          create a group generated by S, with relations L
+createQuotientGroup(L);    create a group generated by the unit matrix whith relations L
+createTorsionFreeGroup(S); create a group generated by S which is torsionfree 
+printGroup(G);             print a group
+
+areIsomorphicGroups(G,H);     test wheter G an H are isomorphic groups (TODO Tuebingen)
+isGroup(G);                   test whether G is a valid group
+isGroupHomomorphism(L1,L2,A); test wheter A defines a group homomrphism from L1 to L2
+
+isGradedRingHomomorphism(R,f,A);  test graded ring homomorph
+createGradedRingHomomorphism(R,f,A);  create a graded ring homomorph
 
 setModuleGrading(M,v);    attach multiweights of units to a module and return it
 getModuleGrading(M);      get multiweights of module units (attached to M)
 
+isSublattice(A,B);        test whether A is a sublattice of B
+imageLattice(P,L);        computes an integral basis for the image of the lattice L under the linear map P.
+intRank(A);               computes the rank of the intmat A
+kernelLattice(P);         computes an integral basis for the kernel of the linear map P.
+latticeBasis(B);          compute an integral basis of the lattice B
+preimageLattice(P,L);     computes an integral basis for the preimage of the lattice L under the linear map P.
+projectLattices(B);       compute a linear map of lattices having the primitive span of B as its kernel.
+intersectLattices(A,B);   compute an integral basis for the intersection of the lattices A and B.
+isIntegralSurjective(P);  test whether the map P of lattices is surjective.
+isPrimitiveSublattice(A); test whether A generates a primitive sublattice. 
+intInverse(A);            compute the integral inverse matrix of the intmat A
+intAdjoint(A,i,j);        delete row i and column j of the intmat A.
+integralSection(P);       for a given linear surjective map P of lattices this procedure returns an integral section of P.
+primitiveSpan(A);         compute a basis for the minimal primitive sublattice that contains the given vectors (by A).
+
+factorgroup(G,H);         create the group G mod H
+productgroup(G,H);        create the group G x H
+
 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
+mDegPartition(p);         compute the multigraded-homogeneous components of p
+
+isTorsionFree();          test whether the current multigrading is  free
+isPositive();             test whether the current multigrading is positive
+isZeroElement(p);         test whether p has zero multidegree
+areZeroElements(M);       test whether an integer matrix M considered as a collection of columns has zero multidegree   
+isHomogeneous(a);         test whether 'a' is multigraded-homogeneous
 
 equalMDeg(e1,e2[,V]);     test whether e1==e2 in the current multigrading
@@ -38 +77 @@
 mDegGroebner(M);          compute the multigraded GB/SB of M
 mDegSyzygy(M);            compute the multigraded syzygies of M
+mDegModulo(I,J);          compute the multigraded 'modulo' module of I and J
 mDegResolution(M,l[,m]);  compute the multigraded resolution of M
-
-defineHomogenous(p);      get a torsion matrix wrt which p becomes homogenous
+mDegTensor(m,n);          compute the tensor product of multigraded modules m,n 
+mDegTor(i,m,n);           compute the Tor_i(m,n) for multigraded modules m,n
+
+defineHomogeneous(p);     get a grading group wrt which p becomes homogeneous
 pushForward(f);           find the finest grading on the image ring, homogenizing f
-
-hermite(A);               compute the Hermite Normal Form of a matrix
+gradiator(h);             coarsens grading of the ring until h becomes homogeneous
+
+hermiteNormalForm(A);     compute the Hermite Normal Form of a matrix
+smithNormalForm(A,#);     compute matrices D,P,Q with D=P*A*Q and D is the smith normal form of A
 
 hilbertSeries(M);         compute the multigraded Hilbert Series of M
+evalHilbertSeries(h,v);   evaluate hilberts series h by substituting v[i] for t_(i)
+
+lll(A);                   applies LLL(.) of lll.lib which only works for lists on a matrix A
 
            (parameters in square brackets [] are optional)
 
-KEYWORDS:  multigradeding, multidegree, multiweights, multigraded-homogenous
+KEYWORDS:  multigrading, multidegree, multiweights, multigraded-homogeneous, integral linear algebra
 ";
 
@@ -56 +103 @@
 
 LIB "standard.lib"; // for groebner
+LIB "lll.lib"; // for lll_matrix
+LIB "matrix.lib"; // for mDegTor
+
+/******************************************************/
+
+static proc concatintmat(intmat A, intmat B)
+{
+
+ if ( nrows(A) != nrows(B) )
+   {
+     ERROR("matrices A and B have different number of rows.");
+   }
+
+ intmat At = transpose(A);
+ intmat Bt = transpose(B);
+
+ intmat Ct[nrows(At) + nrows(Bt)][ncols(At)] = At, Bt;
+
+ return(transpose(Ct));
+}
+
+
+/******************************************************/
+proc createGradedRingHomomorphism(def src, ideal Im, def A)
+"USAGE: createGradedRingHomomorphism(R, f, A); ring R, ideal f, group homomorphism A
+PURPOSE: create a multigraded group ring homomorphism defined by 
+a ring map from R to the current ring, given by generators images f
+and a group homomorphism A between grading groups
+RETURN: graded ring homorphism
+EXAMPLE: example createGradedRingHomomorphism; shows an example
+"
+{
+  string isGRH = "isGRH";
+
+  if( !isGradedRingHomomorphism(def src, ideal Im, def A) ) 
+  {
+    ERROR("Input data is wrong");
+  }
+  
+  list h;
+  h[3] = A;
+  
+//  map f = src, Im;
+  h[2] = Im; // f?
+  h[1] = src;
+  
+  attrib(h, isGRH, (1==1)); // mark it "a graded ring homomorphism"
+
+  return(h);
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+
+  // TODO!
+
+}
+
+
+/******************************************************/
+proc isGradedRingHomomorphism(def src, ideal Im, def A)
+"USAGE: isGradedRingHomomorphism(R, f, A); ring R, ideal f, group homomorphism A
+PURPOSE: test a multigraded group ring homomorphism defined by 
+a ring map from R to the current ring, given by generators images f
+and a group homomorphism A between grading groups
+RETURN: int, 1 for TRUE, 0 otherwise
+EXAMPLE: example isGradedRingHomomorphism; shows an example
+"
+{
+  def dst = basering;
+
+  intmat result_degs = mDeg(Im);
+  print(result_degs);
+   
+  setring src;
+   
+  intmat input_degs = mDeg(maxideal(1));
+  print(input_degs);
+
+  def image_degs = A * input_degs;
+  print( image_degs );
+  
+  def df = image_degs - result_degs;
+  print(df);
+  
+  setring dst;
+   
+  return (areZeroElements( df ));
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+  
+  ring r = 0, (x, y, z), dp;
+  intmat S1[3][3] =
+    1, 0, 0,
+    0, 1, 0,
+    0, 0, 1;
+  intmat L1[3][1] =
+    0,
+    0, 
+    0;
+   
+  def G1 = createGroup(S1, L1); // (S1 + L1)/L1
+  printGroup(G1);
+   
+  setBaseMultigrading(S1, L1); // to change...
+
+  ring R = 0, (a, b, c), dp;
+  intmat S2[2][3] =
+    1, 0,
+    0, 1;
+  intmat L2[2][1] =
+    0,
+    2;
+
+  def G2 = createGroup(S2, L2);
+  printGroup(G2);
+   
+  setBaseMultigrading(S2, L2); // to change...
+
+
+  map F = r, a, b, c;
+  intmat A[nrows(L2)][nrows(L1)] =
+      1, 0, 0,
+      3, 2, -6;
+   
+  // graded ring homomorphism is given by (compatible):
+  print(F);
+  print(A);
+
+  isGradedRingHomomorphism(r, ideal(F), A);
+  def h = createGradedRingHomomorphism(r, ideal(F), A);
+
+  print(h);
+  
+  // not a homo..
+  intmat B[nrows(L2)][nrows(L1)] =
+     1, 1, 1,
+     0, 0, 0;
+  print(B);
+
+  isGradedRingHomomorphism(r, ideal(F), B);
+  def hh = createGradedRingHomomorphism(r, ideal(F), B);
+  
+  if( defined(hh) ) { ERROR("That input was not valid"); }
+}
+
+
+proc createQuotientGroup(intmat L)
+"
+L - relations 
+TODO: bad name 
+"
+{
+  int r = nrows(L); int i;
+  intmat S[r][r]; // SQUARE!!!
+  for(i = r; i > 0; i--){ S[i, i] = 1; }
+  return (createGroup(S,L));
+}
+
+proc createTorsionFreeGroup(intmat S)
+"
+S - generators
+TODO: bad name 
+"
+{
+  int r = nrows(S); int i;
+  intmat L[r][1] = 0;
+  return (createGroup(S,L));
+}
+
+
+/******************************************************/
+proc createGroup(intmat S, intmat L)
+"USAGE: createGroup(S, L); S, L are integer matrices
+PURPOSE: create the group of the form (S+L)/L, i.e. 
+S specifies generators, L specifies relations.
+RETURN: group
+EXAMPLE: example createGroup; shows an example
+"
+{
+  string isGroup = "isGroup";
+  string attrGroupHNF = "hermite";
+  string attrGroupSNF = "smith";
+
+   
+/*
+  if( size(#) > 0 )
+  {
+    if( typeof(#[1]) == "intmat" )
+    {
+      intmat S = #[1];
+    } else { ERROR("Wrong optional argument: 1"); }
+
+    if( size(#) > 1 )
+    {
+      if( typeof(#[2]) == "intmat" )
+      {
+        intmat L = #[2];
+      } else { ERROR("Wrong optional argument: 2"); }
+    }
+  } 
+  if( !defined(S) ) 
+  {}
+*/    
+
+  if( nrows(L) != nrows(S) )
+  {
+    ERROR("Incompatible matrices!");
+  }
+  
+  def H = attrib(L, attrGroupHNF);
+  if( !defined(H) || typeof(H) != "intmat")
+  {
+    attrib(L, attrGroupHNF, hermiteNormalForm(L));
+  } else { kill H; }
+  
+  def HH = attrib(L, attrGroupSNF);
+  if( !defined(HH) || typeof(HH) != "intmat")
+  {
+    attrib(L, attrGroupSNF, smithNormalForm(L));
+  } else { kill HH; }
+
+  list G; // Please, note the order: Generators + Relations:
+  G[1] = S; G[2] = L;
+
+  attrib(G, isGroup, (1==1)); // mark it "a group"
+
+  return (G);
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+   
+  intmat S[3][3] =
+    1, 0, 0,
+    0, 1, 0,
+    0, 0, 1;
+
+  intmat L[3][2] =
+    1, 1,
+    1, 3, 
+    1, 5;
+   
+  def G = createGroup(S, L); // (S+L)/L
+
+  printGroup(G);
+   
+  kill S, L, G;
+
+  /////////////////////////////////////////////////
+  intmat S[2][3] =
+    1, -2, 1,
+    1,  1, 0;
+
+  intmat L[2][1] =
+    0,
+    2;
+
+  def G = createGroup(S, L); // (S+L)/L
+
+  printGroup(G);
+   
+  kill S, L, G;
+   
+  // ----------- extreme case ------------ //
+  intmat S[1][3] =
+    1,  -1, 10;
+
+  // Torsion:
+  intmat L[1][1] =
+    0;
+
+  def G = createGroup(S, L); // (S+L)/L
+
+  printGroup(G);
+}
+
+
+/******************************************************/
+proc printGroup(def G)
+"USAGE: printGroup(G); G is a group
+PURPOSE: prints the group G
+RETURN: nothing
+EXAMPLE: example printGroup; shows an example
+"
+{
+  "Generators: ";
+  print(G[1]);
+
+  "Relations: ";
+  print(G[2]);
+  
+//  attrib(G[2]);
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+   
+}
+
+/******************************************************/
+proc areIsomorphicGroups(def G, def H)
+"USAGE: areIsomorphicGroups(G, H); G and H are groups
+PURPOSE: ?
+RETURN: int, 1 for TRUE, 0 otherwise
+EXAMPLE: example areIsomorphicGroups; shows an example
+"
+{
+  return (1); // TRUE
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+   
+}
+
+
+proc isGroup(def G)
+"test whether G is a valid group"
+{
+  string isGroup = "isGroup";
+  
+  // valid?
+  if( typeof(G) != "list" ){ return(0); }
+
+  def a = attrib(G, isGroup);
+
+///// TODO for Hans: fix attr^2 bug in Singular!
+
+//  if( !defined(a) ) { return(0); }
+//  if( typeof(a) != "int" ) { return(0); } 
+//  if( !a ){ return(0); }
+
+   
+  if( size(G) != 2 ){ return(0); }
+  if( typeof(G[1]) != "intmat" ){ return(0); }
+  if( typeof(G[2]) != "intmat" ){ return(0); }
+  if( nrows(G[1]) != nrows(G[2]) ){ return(0); }
+   
+  return(1==1);
+}
+
+
 
 /******************************************************/
 proc setBaseMultigrading(intmat M, list #)
-"USAGE: setBaseMultigrading(M[, T]); M, T are integer matrices
-PURPOSE: attaches weights of variables and torsion to the basering.
+"USAGE: setBaseMultigrading(M[, G]); M is an intege matrix, G is a group (or lattice)
+PURPOSE: attaches weights of variables and grading group 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
@@ -69 +459 @@
 {
   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);
-  }
+  string attrGradingGroup = "gradingGroup";
+  
+  if( size(#) > 0 )
+  {
+    if( typeof(#[1]) == "intmat" )
+    {
+      def L = createGroup(M, #[1]);
+    } 
+
+    if( isGroup(#[1]) )
+    {
+      def L = #[1];
+
+      if( !isSublattice(M, L[1]) )
+      {
+        ERROR("Multigrading is not contained in the grading group!");
+      }
+    }
+  } 
   else
   {
-    ERROR("Incompatible matrix sizes!");
-  }
+    def L = createTorsionFreeGroup(M);
+  }
+
+  if( !defined(L) ){ ERROR("Wrong arguments: no group given?"); }
+   
+
+  
+  attrib(basering, attrMgrad, M); 
+  attrib(basering, attrGradingGroup, L); 
+
+  ideal Q = ideal(basering);
+  if( !isHomogeneous(Q) ) // easy now, but would be hard before setting ring attributes!
+  {
+    "Warning: your quotient ideal is not homogenous (multigrading was set anyway)!";
+  }
+ 
+   
+
+
 }
 example
 {
   "EXAMPLE:"; echo=2;
-
+   
   ring R = 0, (x, y, z), dp;
 
@@ -106 +512 @@
     0, 0, 1;
 
-  // Torsion:
+  // GradingGroup:
   intmat L[3][2] =
     1, 1,
-    1, 3,
+    1, 3, 
     1, 5;
-
+   
   // attaches M & L to R (==basering):
   setBaseMultigrading(M, L); // Grading: Z^3/L
@@ -117 +523 @@
   // Weights are accessible via "getVariableWeights()":
   getVariableWeights();
-
-  // Test all possible usages:
+  
+  // 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);
+  // Grading group is accessible via "getLattice()":
+  getLattice();
+  
+  // Test all possible usages: 
+  (getLattice() == L) && (getLattice(R) == L) && (getLattice(basering) == L);
+
+  // And its hermite NF via getLattice("hermite"):
+  getLattice("hermite");
+
+  // Test all possible usages: 
+  intmat H = hermiteNormalForm(L);
+  (getLattice("hermite") == H) && (getLattice(R, "hermite") == H) && (getLattice(basering, "hermite") == H);
 
   kill L, M;
@@ -147 +553 @@
     0,
     2;
-
+   
   // attaches M & L to R (==basering):
   setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
@@ -154 +560 @@
   getVariableWeights() == M;
 
-  // Torsion is accessible via "getTorsion()":
-  getTorsion() == L;
+  // Torsion is accessible via "getLattice()":
+  getLattice() == L;
 
   kill L, M;
@@ -167 +573 @@
   intmat L[1][1] =
     0;
-
+   
   // attaches M & L to R (==basering):
   setBaseMultigrading(M); // Grading: Z^3
@@ -174 +580 @@
   getVariableWeights() == M;
 
-  // Torsion is accessible via "getTorsion()":
-  getTorsion() == L;
+  // Torsion is accessible via "getLattice()":
+  getLattice() == L;
 }
 
@@ -208 +614 @@
   def M = attrib(R, attrMgrad);
   if( typeof(M) == "intmat"){ return (M); }
-  ERROR( "Sorry no multigrading matrix!" );
+  ERROR( "Sorry no multigrading matrix!" );  
 }
 example
 {
   "EXAMPLE:"; echo=2;
-
+   
   ring R = 0, (x, y, z), dp;
 
@@ -222 +628 @@
     0, 0, 1;
 
-  // Torsion:
+  // Grading group:
   intmat L[3][2] =
     1, 1,
-    1, 3,
+    1, 3, 
     1, 5;
-
+   
   // attaches M & L to R (==basering):
   setBaseMultigrading(M, L); // Grading: Z^3/L
@@ -243 +649 @@
     1,  1, 0;
 
-  // Torsion:
+  // Grading group:
   intmat L[2][1] =
     0,
     2;
-
+   
   // attaches M & L to R (==basering):
   setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
@@ -262 +668 @@
     1,  -1, 10;
 
-  // Torsion:
+  // Grading group:
   intmat L[1][1] =
     0;
-
+   
   // attaches M & L to R (==basering):
   setBaseMultigrading(M); // Grading: Z^3
@@ -273 +679 @@
 }
 
-/******************************************************/
-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";
+
+proc getGradingGroup(list #)
+"USAGE: getGradingGroup([R])
+PURPOSE: get associated grading group
+RETURN: group, the grading group
+EXAMPLE: example getGradingGroup; shows an example
+"
+{
+  string attrGradingGroup    = "gradingGroup";
 
   int i = 1;
@@ -293 +697 @@
       def R = #[i];
       i++;
-    }
-  }
+    } 
+  } 
 
   if( !defined(R) )
@@ -301 +705 @@
   }
 
-  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);
+  def G = attrib(R, attrGradingGroup);
+  
+  if( !isGroup(G) )
+  {
+    ERROR("Sorry no grading group!");
+  }   
+
+  return(G); 
 }
 example
 {
   "EXAMPLE:"; echo=2;
-
+   
   ring R = 0, (x, y, z), dp;
 
@@ -341 +729 @@
   intmat L[3][2] =
     1, 1,
-    1, 3,
+    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;
+  def G = getGradingGroup();
+
+  printGroup( G );
+  
+  G[1] == M; G[2] == L;
+
+  kill L, M, G;
 
   // ----------- isomorphic multigrading -------- //
@@ -366 +754 @@
     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;
-
+  def G = getGradingGroup();
+
+  printGroup( G );
+  
+  G[1] == M; G[2] == L;
+
+  kill L, M, G;
   // ----------- extreme case ------------ //
 
@@ -387 +774 @@
   intmat L[1][1] =
     0;
-
+   
   // attaches M & L to R (==basering):
   setBaseMultigrading(M); // Grading: Z^3
 
-  // Torsion is accessible via "getTorsion()":
-  getTorsion() == L;
+  def G = getGradingGroup();
+
+  printGroup( G );
+  
+  G[1] == M; G[2] == L;
+
+  kill L, M, G;
+}
+
+
+/******************************************************/
+proc getLattice(list #)
+"USAGE: getLattice([R[,opt]])
+PURPOSE: get associated grading group matrix, i.e. generators (cols) of the grading group
+RETURN: intmat, the grading group matrix, or 
+its hermite normal form if an optional argument (\"hermiteNormalForm\") is given or
+smith normal form if an optional argument (\"smith\") is given
+EXAMPLE: example getLattice; shows an example
+"
+{
+  int i = 1;
+  if( size(#) >= i )
+  {
+    if( ( typeof(#[i]) == "ring" ) or ( typeof(#[i]) == "qring" ) )
+    {
+      i++;
+    } 
+  } 
+
+  string attrGradingGroupHNF = "hermite";
+  string attrGradingGroupSNF = "smith";
+
+  def G = getGradingGroup(#);
+  
+//  printGroup(G); 
+
+
+
+  def T = G[2];
+
+  if( size(#) >= i )
+  {
+    if( #[i] == "hermite" )
+    {
+      def M = attrib(T, attrGradingGroupHNF);
+      if( (!defined(M)) or (typeof(M) != "intmat") )
+      { 
+        M = hermiteNormalForm(T); 
+      }
+      return (M);
+    }
+
+    if( #[i] == "smith" )
+    {
+      def M = attrib(T, attrGradingGroupSNF);
+      if( (!defined(M)) or (typeof(M) != "intmat") )
+      { 
+        M = smithNormalForm(T); 
+      }
+      return (M);
+    }
+  }
+
+  return(T); 
+}
+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 "getLattice()":
+  getLattice() == L;
 
   // its hermite NF:
-  print(getTorsion("hermite"));
-}
+  print(getLattice("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 "getLattice()":
+  getLattice() == L;
+
+  // its hermite NF:
+  print(getLattice("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 "getLattice()":
+  getLattice() == L;
+
+  // its hermite NF:
+  print(getLattice("hermite"));
+}
+
+proc getGradedGenerator(def m, int i)
+"
+returns m[i], but with grading
+"
+{
+  if( typeof(m) == "ideal" ) 
+  {
+    return (m[i]);
+  }
+
+  if( typeof(m) == "module" )
+  {
+    def v = getModuleGrading(m);
+    
+    return ( setModuleGrading(m[i],v) );
+  }
+
+  ERROR("m is expected to be an ideal or a module");
+}
+
 
 /******************************************************/
@@ -410 +948 @@
 
   def V = attrib(m, attrModuleGrading);
-
+  
   if( typeof(V) != "intmat" )
   {
@@ -419 +957 @@
       return (VV);
     }
-
+     
     ERROR("Sorry: vector or module need module-grading-matrix! See 'getModuleGrading'.");
   }
@@ -432 +970 @@
     ERROR("Sorry wrong width of V: " + string(ncols(V)));
   }
-
+   
   return (V);
 }
@@ -438 +976 @@
 {
   "EXAMPLE:"; echo=2;
-
+   
    ring R = 0, (x,y), dp;
    intmat M[2][2]=
@@ -446 +984 @@
      1,  2,  3,  4, 0,
      0, 10, 20, 30, 1;
-
+   
    setBaseMultigrading(M, T);
-
+   
    ideal I = x, y, xy^5;
-   isHomogenous(I);
-
+   isHomogeneous(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);
+   
+   vector v = getGradedGenerator(S, 1);
    getModuleGrading(v) == V;
-   isHomogenous(v);
-   print( mDeg(v) );
-
-   isHomogenous(S);
+   isHomogeneous(v);
+   print( mDeg(v) );   
+   
+   isHomogeneous(S);
    print( mDeg(S) );
 }
@@ -474 +1012 @@
 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).
+         module/vector homogeneous. One can do that using isHomogeneous(m).
 EXAMPLE: example setModuleGrading; shows an example
 "
@@ -484 +1022 @@
   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);
@@ -491 +1029 @@
 {
   "EXAMPLE:"; echo=2;
-
+   
    ring R = 0, (x,y), dp;
    intmat M[2][2]=
@@ -499 +1037 @@
      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);
+   
+   vector v = getGradedGenerator(S, 1);
    getModuleGrading(v) == V;
 
-   print( mDeg(v) );
-
-   isHomogenous(S);
+   print( mDeg(v) );   
+
+   isHomogeneous(S);
 
    print( mDeg(S) );
@@ -526 +1064 @@
 
 
-
+proc mDegTensor(module m, module n){
+  matrix M = m;
+  matrix N = n;
+  intmat gm = getModuleGrading(m);
+  intmat gn = getModuleGrading(n);
+  int grows = nrows(gm);
+  int mr = nrows(M);
+  int mc = ncols(M);
+  if(rank(M) == 0){ mc = 0;}
+  int nr = nrows(N);
+  int nc = ncols(N);
+  if(rank(N) == 0){ nc = 0;}
+  intmat gresult[nrows(gm)][mr*nr];
+  matrix result[mr*nr][mr*nc+mc*nr];
+  int i, j;
+  int column = 1;
+  for(i = 1; i<=mr; i++){
+    for(j = 1; j<=nr; j++){
+      gresult[1..grows,(i-1)*nr+j] = gm[1..grows,i]+gn[1..grows,j];
+    }
+  }
+  //gresult;
+  if( nc!=0 ){
+    for(i = 1; i<=mr; i++)
+    {
+      result[((i-1)*nr+1)..(i*nr),((i-1)*nc+1)..(i*nc)] = N[1..nr,1..nc]; 
+    }
+  }
+  list rownumbers, colnumbers;
+  //print(result);
+  if( mc!=0 ){
+    for(j = 1; j<=nr; j++)
+    {
+      rownumbers = nr*(0..(mr-1))+j*(1:mr);
+      colnumbers = ((mr*nc+j):mc)+nr*(0..(mc-1));
+      result[rownumbers[1..mr],colnumbers[1..mc] ] = M[1..mr,1..mc];
+    }
+  }
+  module res = result;
+  res = setModuleGrading(res, gresult);
+  //getModuleGrading(res);
+  return(res);
+}
+example
+{
+"EXAMPLE: ";echo=2;
+ring r = 0,(x),dp;
+intmat g[2][1]=1,1;
+setBaseMultigrading(g);
+matrix m[5][3]=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15;
+matrix n[3][2]=x,x2,x3,x4,x5,x6;
+module mm = m;
+module nn = n;
+intmat gm[2][5]=1,2,3,4,5,0,0,0,0,0;
+intmat gn[2][3]=0,0,0,1,2,3;
+mm = setModuleGrading(mm, gm);
+nn = setModuleGrading(nn, gn);
+module mmtnn = mDegTensor(mm, nn);
+print(mmtnn);
+getModuleGrading(mmtnn);
+LIB "homolog.lib";
+module tt = tensorMod(mm,nn);
+print(tt);
+
+kill m, mm, n, nn, gm, gn;
+
+matrix m[7][3] = x, x-1,x+2, 3x, 4x, x5, x6, x-7, x-8, 9, 10, 11x, 12 -x, 13x, 14x, x15, (x-4)^2, x17, 18x, 19x, 20x, 21x;
+matrix n[2][4] = 1, 2, 3, 4, x, x2, x3, x4;
+module mm = m;
+module nn = n;
+print(mm);
+print(nn);
+intmat gm[2][7] = 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 0;
+intmat gn[2][2] = 0, 0, 1, 2;
+mm = setModuleGrading(mm, gm);
+nn = setModuleGrading(nn, gn);
+module mmtnn = mDegTensor(mm, nn);
+print(mmtnn);
+getModuleGrading(mmtnn);
+matrix a = mmtnn;
+matrix b = tensorMod(mm, nn);
+print(a-b);
+
+}
+
+proc mDegTor(int i, module m, module n)
+{
+  def res = mDegResolution(n, 0, 1);
+  //print(res);
+  list l = res;
+  if(size(l)<i){ return(0);}
+  else
+  {
+    
+    matrix fd[nrows(m)][0];
+    matrix fd2[nrows(l[i+1])][0];
+    matrix fd3[nrows(l[i])][0];
+
+    module freedim = fd;
+    module freedim2 = fd2;
+    module freedim3 = fd3;
+
+    freedim = setModuleGrading(freedim,getModuleGrading(m));
+    freedim2 = setModuleGrading(freedim2,getModuleGrading(l[i+1]));
+    freedim3 = setModuleGrading(freedim3, getModuleGrading(l[i]));
+    
+    module mimag = mDegTensor(freedim3, m);
+    //"mimag ok.";
+    module mf = mDegTensor(l[i], freedim);
+    //"mf ok.";
+    module mim1 = mDegTensor(freedim2 ,m);
+    module mim2 = mDegTensor(l[i+1],freedim);
+    //"mim1+2 ok.";
+    module mker = mDegModulo(mf,mimag);
+    //"mker ok.";
+    module mim = mim1,mim2;
+    mim = setModuleGrading(mim, getModuleGrading(mim1));
+    //"mim: r: ",nrows(mim)," c: ",ncols(mim);
+    //"mim1: r: ",nrows(mim1)," c: ",ncols(mim1);
+    //"mim2: r: ",nrows(mim2)," c: ",ncols(mim2);
+    //matrix mimmat = mim;
+    //matrix mimmat1[16][4]=mimmat[1..16,25..28];
+    //print(mimmat1-matrix(mim2));
+    return(mDegModulo(mker,mim));
+    //return(0);
+  }
+  return(0);
+}
+example
+{
+"EXAMPLE: ";echo=2;
+LIB "homolog.lib";
+ring r = 0,(x_(1..4)),dp;
+intmat g[2][4]=1,1,0,0,0,1,1,-1;
+setBaseMultigrading(g);
+ideal i = maxideal(1);
+module m = mDegSyzygy(i);
+module rt = Tor(2,m,m);
+module mDegT = mDegTor(2,m,m);
+print(matrix(rt)-matrix(mDegT));
+/*
+ring r = 0,(x),dp;
+intmat g[2][1]=1,1;
+setBaseMultigrading(g);
+matrix m[5][3]=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15;
+matrix n[3][2]=x,x2,x3,x4,x5,x6;
+module mm = m;
+module nn = n;
+intmat gm[2][5]=1,1,1,1,1,1,1,1,1,1,1;
+intmat gn[2][3]=0,-2,-4,0,-2,-4;
+mm = setModuleGrading(mm, gm);
+nn = setModuleGrading(nn, gn);
+isHomogeneous(mm,"checkGens");
+isHomogeneous(nn,"checkGens");
+mDegTor(1,mm, nn);
+
+kill m, mm, n, nn, gm, gn;
+
+matrix m[7][3] = x, x-1,x+2, 3x, 4x, x5, x6, x-7, x-8, 9, 10, 11x, 12 -x, 13x, 14x, x15, (x-4)^2, x17, 18x, 19x, 20x, 21x;
+matrix n[2][4] = 1, 2, 3, 4, x, x2, x3, x4;
+module mm = m;
+module nn = n;
+print(mm);
+print(nn);
+intmat gm[2][7] = 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 0;
+intmat gn[2][2] = 0, 0, 1, 2;
+mm = setModuleGrading(mm, gm);
+nn = setModuleGrading(nn, gn);
+module mmtnn = mDegTensor(mm, nn);
+*/
+}
+
+
+/******************************************************/
+proc isGroupHomomorphism(def L1, def L2, intmat A)
+"USAGE: gisGoupHomomorphism(L1,L2,A); L1 and L2 are groups, A is an integer matrix
+PURPOSE: checks whether A defines a group homomorphism phi: L1 --> L2
+RETURN: int, 1 if A defines the homomorphism and 0 otherwise 
+EXAMPLE: example isGroupHomomorphism; shows an example
+"
+{
+  // TODO: L1, L2 
+  if( (ncols(A) != nrows(L1)) or (nrows(A) != nrows(L2)) )
+  {
+    ERROR("Incompatible sizes!");
+  }
+
+  intmat im = A * L1;
+  
+  return  (areZeroElements(im, L2)); 
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+
+   intmat L1[4][1]=
+     0,
+     0,
+     0,
+     2;
+     
+   intmat L2[3][2]=
+     0, 0,
+     2, 0,
+     0, 3;
+
+  intmat A[3][4] = 
+     1, 2, 3, 0,
+     7, 0, 0, 0,
+     1, 2, 0, 3;
+  print( A );
+
+  isGroupHomomorphism(L1, L2, A);
+
+  intmat B[3][4] = 
+     1, 2, 3, 0,
+     7, 0, 0, 0,
+     1, 2, 0, 2;
+  print( B ); 
+
+  isGroupHomomorphism(L1, L2, B); // Not a homomorphism!
+}
 
 /******************************************************/
 proc isTorsionFree()
 "USAGE: isTorsionFree()
-PURPOSE: Determines whether the multigrading attached to the current ring is torsion-free.
+PURPOSE: Determines whether the multigrading attached to the current ring is free.
 RETURN: boolean, the result of the test
 EXAMPLE: example isTorsionFree; shows an example
 "
 {
-
-  intmat H = hermite(transpose(getTorsion("hermite"))); // TODO: ?cache it?  //******
+  intmat H = smithNormalForm(getLattice()); // TODO: ?cache it?  //******
 
   int i, j;
@@ -550 +1308 @@
 
     if(H[j, i]!=0)
-    {
+    { 
       d=d*H[j, i];
     }
@@ -556 +1314 @@
 
   if( (d*d)==1 )
-  {
+  { 
     return(1==1);
   }
@@ -572 +1330 @@
      1, 2, 3, 4, 0,
      0,10,20,30, 1;
-
+   
    setBaseMultigrading(M,T);
-
-   // Is the resulting group torsion free?
+   
+   // Is the resulting group  free?
    isTorsionFree();
 
@@ -582 +1340 @@
 
    ring R=0,(x,y,z),dp;
-   intmat A[3][3] =
+   intmat A[3][3] = 
      1,0,0,
      0,1,0,
@@ -591 +1349 @@
      1,2,0,3;
    setBaseMultigrading(A,B);
-   // Is the resulting group torsion free?
+   // Is the resulting group  free?
    isTorsionFree();
 
@@ -598 +1356 @@
 
 
+static proc gcdcomb(int a, int b)
+{
+  // a;
+  // b;
+  intvec av = a,1,0;
+  intvec bv = b,0,1;
+  intvec save;
+  while(av[1]*bv[1] != 0)
+  {
+    bv = bv - (bv[1] - bv[1]%av[1])/av[1] * av;
+    save = bv;
+    bv = av;
+    av = save;
+  }
+  if(bv[1] < 0)
+  {
+    bv = -bv;
+  }
+  return(bv);
+}
+
+
+proc lll(def A)
+"
+The lll algorithm of lll.lib only works for lists of vectors.
+Maybe one should rescript it for matrices. This method will
+convert a matrix to a list, plug it into lll and make the result
+a matrix and return it.
+"
+{
+  if(typeof(A) == "list")
+  {
+    int sizeA= size (A);
+    if (sizeA == 0)
+    {
+      return (A);
+    }
+    if (typeof (A [1]) != "intvec")
+    {
+      ERROR("Unrecognized type.");
+    }
+    int columns= size (A [1]);
+    int i;
+    for (i= 2; i <= sizeA; i++)
+    {
+      if (typeof (A[i]) != "intvec")
+      {
+        ERROR("Unrecognized type.");
+      }
+      if (size (A [i]) != columns)
+      {
+        ERROR ("expected equal dimension");
+      }
+    }
+    int j;
+    intmat m [columns] [sizeA];
+    for (i= 1; i <= sizeA; i++)
+    {
+      for (j= 1; j <= columns; j++)
+      {
+        m[i,j]= A[i] [j];
+      }
+    }
+    m= system ("LLL", m);
+    list result= list();
+
+    for (i= 1; i <= sizeA; i++)
+    {
+      intvec buf= intvec (m[i , 1]);
+      for (j= 2; j <= columns; j++)
+      {
+        buf= buf,intvec (m [i, j]);
+      }
+      result= result+ list (buf);
+      
+    }
+    return(result);
+  } 
+  else
+  {
+    if(typeof(A) == "intmat")
+    {
+      A= system ("LLL", A);
+      return(A);
+    }
+    else
+    {
+      ERROR("Unrecognized type.");
+    }
+  }
+}
+
+example
+{
+
+"EXAMPLE:";
+
+ring R = 0,x,dp;
+intmat m[5][5]=13,25,37,83,294,12,-33,9,0,64,77,12,34,6,1,43,2,88,91,100,-46,32,37,42,15;
+lll(m);
+list l=intvec(13,25,37, 83, 294),intvec(12, -33, 9,0,64), intvec (77,12,34,6,1), intvec (43,2,88,91,100), intvec (-46,32,37,42,15);
+lll(l);
+}
+
+
+proc smithNormalForm(intmat A, list #)
+"
+This method returns 3 Matrices P, D and Q such that D = P*A*Q.
+WARNING: This might not be what you expect.
+"
+{
+  list l1 = hermiteNormalForm(A, 5);
+  // l1;
+  intmat B = transpose(l1[1]);
+  list l2 = hermiteNormalForm(B, 5);
+  // l2;
+  intmat P = transpose(l2[2]);
+  intmat D = transpose(l2[1]);
+  intmat Q = l1[2];
+  int cc = ncols(D);
+  int rr = nrows(D);
+  intmat transform;
+  int k = 1;
+  int a, b, c;
+  // D;
+  intvec v;
+  if((cc==1)||(rr==1)){
+    if(size(#)==0)
+    {
+      return(D);
+    } else
+    {
+      return(list(P,D,Q));
+    }
+  }
+  while(D[k+1,k+1] !=0){
+    if(D[k+1,k+1]%D[k,k]!=0){
+      b = D[k, k]; c = D[k+1, k+1];
+      v = gcdcomb(D[k,k],D[k+1,k+1]);
+      transform = unitMatrix(cc);
+      transform[k+1,k] = 1;
+      a = -v[3]*D[k+1,k+1]/v[1];
+      transform[k, k+1] = a;
+      transform[k+1, k+1] = a+1;
+      //det(transform);
+      D = D*transform;
+      Q = Q*transform;
+      //D;
+      transform = unitMatrix(rr);
+      transform[k,k] = v[2];
+      transform[k,k+1] = v[3];
+      transform[k+1,k] = -c/v[1];
+      transform[k+1,k+1] = b/v[1];
+      D = transform * D;
+      P = transform * P;
+      //" ";
+      //D;
+      //"small transform: ", det(transform);
+      //transform;
+      k=0;
+    }
+    k++;
+    if((k==rr) || (k==cc)){
+      break;
+    }
+  }
+  //"here is the size ",size(#); 
+  if(size(#) == 0){
+    return(D);
+  } else {
+    return(list(P, D, Q));
+  }
+}
+example
+{
+"EXAMPLE: "; echo=2;
+
+intmat A[5][7] = 
+1,0,1,0,-2,9,-71,
+0,-24,248,-32,-96,448,-3496,
+0,4,-42,4,-8,30,-260,
+0,0,0,18,-90,408,-3168,
+0,0,0,-32,224,-1008,7872;
+
+list l = smithNormalForm(A, 5);
+
+l;
+l[1]*A*l[3];
+det(l[1]);
+det(l[3]);
+}
+
 
 /******************************************************/
-proc hermite(intmat A)
-"USAGE: hermite( A );
-PURPOSE: Computes the (lower triangular) Hermite Normal Form
+proc hermiteNormalForm(intmat A, list #) 
+"USAGE: hermiteNormalForm( 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)
+EXAMPLE: example hermiteNormalForm; shows an example
+"
+{
+  
+  int row, column, i, j;
+  int rr = nrows(A);
+  int cc = ncols(A);
+  intvec savev, gcdvec, v1, v2;
+  intmat q = unitMatrix(cc);
+  intmat transform;
+  column = 1;
+  for(row = 1; (row<=rr)&&(column<=cc); row++)
+    {
+      if(A[row,column]==0)
         {
-          for(k=1;k<=nrows(A);k++)
-          {
-            save=A[k, i];
-            A[k, i]=A[k, j];
-            A[k, j]=save;
-          }
-          break;
+          for(j = column; j<=cc; j++)
+            {
+              if(A[row, j]!=0)
+                {
+                  transform = unitMatrix(cc);
+                  transform[j,j] = 0;
+                  transform[column, column] = 0;
+                  transform[column,j] = 1;
+                  transform[j,column] = 1;
+                  q = q*transform;
+                  A = A*transform;
+                  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)
+      if(A[row,column] == 0)
         {
-          for(k=1;k<=nrows(A);k++){ A[k, j]=(-1)*A[k, j];}
+          row++;
+          continue;
         }
-
-        if(A[i+l, i]==1)
+      for(j = column+1; j<=cc; j++)
         {
-          d=A[i+l, j];
-          for(k=1;k<=nrows(A);k++)
-          {
-            A[k, j]=A[k, j]-d*A[k, i];
-          }
+          if(A[row, j]!=0)
+            {
+              gcdvec = gcdcomb(A[row,column],A[row,j]);
+              // gcdvec;
+              // typeof(A[1..rr,column]);
+              v1 = A[1..rr,column];
+              v2 = A[1..rr,j];
+              transform = unitMatrix(cc);
+              transform[j,j] = v1[row]/gcdvec[1];
+              transform[column, column] = gcdvec[2];
+              transform[column,j] = -v2[row]/gcdvec[1];
+              transform[j,column] = gcdvec[3];
+              q = q*transform;
+              A = A*transform;
+              // A;
+            }
         }
-        else
+      if(A[row,column]<0)
         {
-          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];
-            }
-          }
+          transform = unitMatrix(cc);
+          transform[column,column] = -1;
+          q = q*transform;
+          A = A*transform;
+        }
+      for( j=1; j<column; j++){
+        if(A[row, j]!=0){
+          transform = unitMatrix(cc);
+          transform[column, j] = (-A[row,j]+A[row, j]%A[row, column])/A[row, column];
+          if(A[row,j]<0){
+            transform[column,j]=transform[column,j]+1;}
+          q = q*transform;
+          A = A*transform; 
         }
       }
-      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);
+      column++;
+    }
+  if(size(#) > 0){
+    return(list(A, q));
+  }
+  return(A);
 }
 example
@@ -697 +1637 @@
   "EXAMPLE:"; echo=2;
 
-   intmat M[2][5] =
+   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] =
+   print(hermiteNormalForm(M));
+
+   intmat T[3][4] = 
      3,3,3,3,
      2,1,3,0,
@@ -710 +1650 @@
 
    // Hermite Normal Form of T:
-   print(hermite(T));
-
-   intmat A[4][5] =
+   print(hermiteNormalForm(T));
+
+   intmat A[4][5] = 
      1,2,3,2,2,
      1,2,3,4,0,
@@ -719 +1659 @@
 
    // 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 );
-
-}
+   print(hermiteNormalForm(A));
+}
+
+proc areZeroElements(intmat m, list #)
+"same as isZeroElement but for an integer matrix considered as a collection of columns"
+{
+  int r = nrows(m);
+  int i = ncols(m);
+
+  intvec v;
+  
+  for( ; i > 0; i-- )
+  {
+    v = m[1..r, i];
+    if( !isZeroElement(v, #) )
+    { 
+      return (0);
+    }
+  }
+  return(1);
+}
+
 example
 {
@@ -776 +1694 @@
     1;
 
+  intmat tt[2][1]=
+    1,
+    -1;
+
   setBaseMultigrading(g,t);
 
@@ -785 +1707 @@
   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);
+ intmat m[5][2]=mDeg(a)-mDeg(b),mDeg(b)-mDeg(c),mDeg(c)-mDeg(d),mDeg(d)-mDeg(e),mDeg(e)-mDeg(f);
+ m=transpose(m);
+ areZeroElementes(m);
+ areZeroElementes(m,tt);
 }
 
 
 /******************************************************/
-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
+proc isZeroElement(intvec mdeg, list #)
+"USAGE: isZeroElement(d, [T]); intvec d, group T
+PURPOSE: For a integer vector mdeg representing the multidegree of some polynomial
+or vector this method computes if the multidegree is contained in the grading group
+group (either set globally or given as an optional argument), i.e. if it is zero in the multigrading.
+EXAMPLE: example isZeroElement; shows an example
 "
 {
@@ -815 +1727 @@
     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);
+      intmat H = hermiteNormalForm(#[1]);
+    } else
+    {
+      if( typeof(#[1]) == "list" )
+      {
+        list L = #[1];
+        intmat H = attrib(L, "hermite"); // todo
+      }
+    }
+    
+  }
+  if( !defined(H) )
+  {
+    intmat H = getLattice("hermite");
+  }
+
+  int x, k, i, row;
+
+  int r = nrows(H);
+  int c = ncols(H);
+
+  int rr = nrows(mdeg);
+  row = 1;
   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;
+  for(i=1; (i<=r)&&(row<=r)&&(i<=c); i++)
+  {
+    while((H[row,i]==0)&&(row<=r))
+    {
+      row++;
+      if(row == (r+1)){
+        break;
+      }
+    }
+    if(row<=r){
+      if(H[row,i]!=0)
+      {
+        v = H[1..r,i];
+        mdeg = mdeg-(mdeg[row]-mdeg[row]%v[row])/v[row]*v;
+      }
+    }
+  }
+  return( mdeg == 0 );
+
+}
+example
+{
+ "EXAMPLE:"; echo=2;
 
   ring r = 0,(x,y,z),dp;
@@ -1075 +1780 @@
     1,0,1,
     0,1,1;
-
   intmat t[2][1]=
     -2,
     1;
+
+  intmat tt[2][1]=
+    1,
+    -1;
 
   setBaseMultigrading(g,t);
@@ -1089 +1797 @@
   poly f = z2;
 
+  intvec v1 = mDeg(a) - mDeg(b);
+  v1;
+  isZeroElement(v1);
+  isZeroElement(v1, tt);
+  
+  intvec v2 = mDeg(a) - mDeg(c);
+  v2;
+  isZeroElement(v2);
+  isZeroElement(v2, tt);
+  
+  intvec v3 = mDeg(e) - mDeg(f);
+  v3;
+  isZeroElement(v3);
+  isZeroElement(v3, tt);
+  
+  intvec v4 = mDeg(c) - mDeg(d);
+  v4;
+  isZeroElement(v4);
+  isZeroElement(v4, tt);
+}
+
+
+/******************************************************/
+proc defineHomogeneous(poly f, list #)
+"USAGE: defineHomogeneous(f[, G]); polynomial f, integer matrix G
+PURPOSE: Yields a matrix which has to be appended to the grading group matrix to make the 
+polynomial f homogeneous in the grading by grad.
+EXAMPLE: example defineHomogeneous; shows an example
+"
+{
+  if( size(#) > 0 ) 
+  {
+    if( typeof(#[1]) == "intmat" ) 
+    {
+      intmat grad = #[1];
+    }
+  }
+
+  if( !defined(grad) )
+  {
+    intmat grad = getVariableWeights();
+  }
+
+  intmat newgg[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++)
+    {
+      newgg[j,i-1] = v[j];
+    }
+  }
+  return(newgg);
+}
+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 = defineHomogeneous(f);
+  M;
+  defineHomogeneous(f, grad) == M;
+  
+  isHomogeneous(f);
+  setBaseMultigrading(grad, M);
+  isHomogeneous(f);
+}
+
+
+proc gradiator(def h)
+PURPOSE: coarsens the grading of the basering until the polynom or ideal h becomes homogeneous.
+
+{
+  if(typeof(h)=="poly"){
+    intmat W = getVariableWeights();
+    intmat L = getLattice();
+    intmat toadd = defineHomogeneous(h);
+    //h;
+    //toadd;
+    if(ncols(toadd) == 0)
+    {
+      return(1==1);
+    }
+    int rr = nrows(W);
+    intmat newL[rr][ncols(L)+ncols(toadd)];
+    newL[1..rr,1..ncols(L)] = L[1..rr,1..ncols(L)];
+    newL[1..rr,(ncols(L)+1)..(ncols(L)+ncols(toadd))] = toadd[1..rr,1..ncols(toadd)];
+    setBaseMultigrading(W,newL);
+    return(1==1);
+  }
+  if(typeof(h)=="ideal"){
+    int i;
+    def s = (1==1);
+    for(i=1;i<=size(h);i++){
+      s = s && gradiator(h[i]);
+    }
+    return(s);
+  }
+  return(1==0); 
+}
+example
+{
+"EXAMPLE:"; echo=2;
+ring r = 0,(x,y,z),dp;
+intmat g[2][3] = 1,0,1,0,1,1;
+intmat l[2][1] = 3,0;
+
+setBaseMultigrading(g,l);
+
+getLattice();
+
+ideal i = -y5+x4,
+          y6+xz,
+          x2y;
+gradiator(i);
+getLattice();
+isHomogeneous(i);
+}
+
+
+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 grading group 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=getLattice(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 = defineHomogeneous(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 = defineHomogeneous(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 gragr[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++){gragr[i,j+columns]=newtor[k][i,j];}
+    }
+    columns=columns+ncols(newtor[k]);
+  }
+
+  gragr=hermiteNormalForm(gragr);
+  intmat result[pp][pp];
+  for(i=1;i<=pp;i++){
+    for(j=1;j<=pp;j++){result[i,j]=gragr[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();
+  getLattice();
+
+  // 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();
+  getLattice();
+  
+
+  // 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 grading group!?
+  {
+    return ( difference == 0);
+  }
+  return ( isZeroElement( difference ) );
+}
+example
+{
+  "EXAMPLE:"; echo=2;printlevel=3;
+
+  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));
@@ -1116 +2183 @@
 /******************************************************/
 static proc isFreeRepresented()
-"check whether the base muligrading is torsion-free (it is zero).
-"
-{
-  intmat T = getTorsion();
+"check whether the base muligrading is free (it is zero).
+"
+{
+  intmat T = getLattice();
 
   intmat Z[nrows(T)][ncols(T)];
 
-  return (T == Z); // no torsion!
+  return (T == Z); // no grading group!
 }
 
 
 /******************************************************/
-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
+proc isHomogeneous(def a, list #)
+"USAGE: isHomogeneous(a[, f]); a polynomial/vector/ideal/module
+RETURN: boolean, TRUE if a is (multi)homogeneous, and FALSE otherwise
+EXAMPLE: example isHomogeneous; shows an example
 "
 {
@@ -1137 +2204 @@
   {
     return ( size(mDegPartition(a)) <= 1 )
-  }
-
-  if( typeof(a) == "vector" or typeof(a) == "module" )
-  {
-    intmat V = getModuleGrading(a);
   }
 
@@ -1153 +2215 @@
         for( int i = ncols(a); i > 0; i-- )
         {
-          aa = a[i];
-
-          if(defined(V))
-          {
-            aa = setModuleGrading(aa, V);
-          }
-
-          if(!isHomogenous(aa))
+          aa = getGradedGenerator(a, i);
+
+          if(!isHomogeneous(aa))
           {
             return(0==1);
@@ -1171 +2228 @@
     def g = groebner(a); // !!!!
 
-    def b, aa; int j;
+    def b, aa; int j; 
     for( int i = ncols(a); i > 0; i-- )
     {
-      aa = a[i];
-
-      if(defined(V))
-      {
-        aa = setModuleGrading(aa, V);
-      }
+      aa = getGradedGenerator(a, i);
 
       b = mDegPartition(aa);
@@ -1191 +2243 @@
     }
     return(1==1);
-  }
+  }  
 }
 example
@@ -1200 +2252 @@
 
   //Grading and Torsion matrices:
-  intmat M[3][3] =
+  intmat M[3][3] = 
     1,0,0,
     0,1,0,
@@ -1217 +2269 @@
   print(mDeg(_));
 
-  isHomogenous(f);   // f: is not homogenous
+  isHomogeneous(f);   // f: is not homogeneous
 
   poly g = 1-xy2z3;
-  isHomogenous(g); // g: is homogenous
+  isHomogeneous(g); // g: is homogeneous
   mDegPartition(g);
 
@@ -1226 +2278 @@
   /////////////////////////////////////////////////////////
   // new Torsion matrix:
-  intmat T[3][4] =
+  intmat T[3][4] = 
     3,3,3,3,
     2,1,3,0,
     1,2,0,3;
-
+  
   setBaseMultigrading(M,T);
 
   f;
-  isHomogenous(f);
+  isHomogeneous(f);
   mDegPartition(f);
 
-  // ---------------------
+  // --------------------- 
   g;
-  isHomogenous(g);
+  isHomogeneous(g);
   mDegPartition(g);
 
@@ -1246 +2298 @@
   ring R = 0, (x,y,z), dp;
 
-  intmat A[2][3] =
+  intmat A[2][3] = 
     0,0,1,
     3,2,1;
-  intmat T[2][1] =
-    -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
+  isHomogeneous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3)); // 1
+  isHomogeneous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3), "checkGens");
+  isHomogeneous(ideal(x+y, x2 - y2)); // 0
 
   // Degree partition:
@@ -1262 +2314 @@
   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);
@@ -1278 +2330 @@
 
   N = setModuleGrading(N, V);
-  isHomogenous(N);
+  isHomogeneous(N);
   print( mDeg(N) );
 
-  ///////////////////////////////////////
-
-  V =
-    1, 2, 3,
+  /////////////////////////////////////// 
+
+  V = 
+    1, 2, 3,  
     4, 5, 6;
 
@@ -1296 +2348 @@
 
   N = setModuleGrading(N, V);
-  isHomogenous(N);
+  isHomogeneous(N);
   print( mDeg(N) );
 
-  ///////////////////////////////////////
-
-  V =
-    0, 0, 0,
+  /////////////////////////////////////// 
+
+  V = 
+    0, 0, 0,  
     4, 1, 0;
 
   N = gen(1) + x * gen(2), z*gen(3);
   N = setModuleGrading(N, V); print(N);
-  isHomogenous(N);
+  isHomogeneous(N);
   print( mDeg(N) );
-  v1 = setModuleGrading(N[1], V); print(v1);
+  v1 = getGradedGenerator(N,1); print(v1);
   mDegPartition(v1);
   print( mDeg(_) );
   N = setModuleGrading(N, V); print(N);
-  isHomogenous(N);
+  isHomogeneous(N);
   print( mDeg(N) );
 }
@@ -1339 +2391 @@
   int r = nrows(M);
 
+  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)));
+    }
+  }
+                            
   if( A == 0 )
   {
@@ -1345 +2407 @@
   }
 
-  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;
 
@@ -1367 +2419 @@
     A = A - lead(A);
     while( size(A) > 0 )
-    {
+    { 
       v = leadexp(A); //  v;
       m = max( m, M * v, r ); // ????
@@ -1388 +2440 @@
     A = A - lead(A);
     while( size(A) > 0 )
-    {
+    { 
       v = leadexp(A); //  v;
 
@@ -1413 +2465 @@
       {
         G[j, i] = d[j];
-      }
+      }      
     }
     return(G);
@@ -1425 +2477 @@
     for( i = ncols(A); i > 0; i-- )
     {
-      v = setModuleGrading(A[i], V);
-
-      // G[1..r, i]
+      v = getGradedGenerator(A, i); 
+
+      // G[1..r, i] 
       d = mDeg(v);
 
@@ -1433 +2485 @@
       {
         G[j, i] = d[j];
-      }
+      }      
 
     }
@@ -1439 +2491 @@
     return(G);
   }
-
+ 
 }
 example
@@ -1453 +2505 @@
   print(Ta);
 
-  //   attrib(A, "torsion", Ta); // to think about
+  //   attrib(A, "gradingGroup", Ta); // to think about
 
 //  "poly:";
@@ -1463 +2515 @@
 
   setBaseMultigrading(A, Ta);
-
+  
   mDeg( x*x*y );
-
+  
   mDeg( y*y*y*x );
-
+  
   mDeg( x*y + x + 1 );
 
@@ -1477 +2529 @@
 
 //  "ideal:";
-
+  
   ideal I = y*x*x, x*y*y*y;
   print( mDeg(I) );
@@ -1485 +2537 @@
 
 //  "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);
@@ -1499 +2551 @@
   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];
@@ -1518 +2570 @@
   DDD = setModuleGrading(DDD, B);
   print( mDeg( DDD ) );
-}
+
+};
+
+
+
+
 
 /******************************************************/
@@ -1526 +2583 @@
 EXAMPLE: example mDegPartition; shows an example
 "
-{
-  if( typeof(p) == "poly" )
-  {
-    ideal I;
+{ // TODO: What about an ideal or module???
+
+  if( typeof(p) == "poly" ) 
+  {
+    ideal I;      
     poly mp, t, tt;
   }
@@ -1536 +2594 @@
     if(  typeof(p) == "vector" )
     {
-      module I;
+      module I;      
       vector mp, t, tt;
     }
@@ -1547 +2605 @@
   if( typeof(p) == "vector" )
   {
-    intmat V = getModuleGrading(p);
+    intmat V = getModuleGrading(p); 
   }
   else
@@ -1554 +2612 @@
   }
 
-  if( size(p) > 1)
+  if( size(p) > 1) 
   {
     intvec m;
@@ -1561 +2619 @@
     {
       m = leadexp(p);
-      mp = lead(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  ) )
+        if( equalMDeg( leadexp(tt), m, V  ) ) 
         {
           mp = mp + lead(tt); // "mp", mp;
@@ -1616 +2674 @@
 
   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);
 }
@@ -1631 +2689 @@
 {
   intmat A[n][n];
-
+  
   for( int i = n; i > 0; i-- )
   {
@@ -1646 +2704 @@
 "
 USAGE: finestMDeg(r); ring r
-RETURN: ring, r endowed with the finest multigrading
+RETURN: ring, r endowed with the finest multigrading 
 TODO: not yet...
 "
@@ -1675 +2733 @@
 
 
-  if( n > 0)
-  {
-
-    intmat L[N][n];
+  if( n > 0) 
+  {
+
+    intmat L[N][n]; 
     //  list L;
     int j = n;
@@ -1686 +2744 @@
       p = I[i];
 
-      if( size(p) > 1 )
+      if( size(p) > 1 ) 
       {
         intvec m0 = leadexp(p);
@@ -1701 +2759 @@
 
     print(L);
-    setBaseMultigrading(A, L);
-  }
+    setBaseMultigrading(A, L);      
+  } 
   else
   {
@@ -1739 +2797 @@
 
   if( size(#) > 0 and typeof(#[1]) == "intmat" )
-  {
+  { 
     attrib(F, "P", #[1]);
   }
@@ -1797 +2855 @@
   // check:
   print(L*M);
-}
+};
+
+
 
 /******************************************************/
@@ -1819 +2879 @@
 "
 {
+   if( system("sh","which hilbert 2> /dev/null 1> /dev/null") != 0 )
+   {
+     ERROR("Sorry: cannot find 'hilbert' command from 4ti2. Please install 4ti2!");
+   }
+   
 //--------------------------------------------------------------------------
 // Initialization and Sanity Checks
@@ -1870 +2935 @@
 // 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 " +
+   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
@@ -1889 +2954 @@
 
    string s = read(ausg);
+
+   if( defined(keepfiles) <= 0)
+   {
+      j=system("sh",("rm -f sing4ti2.hil sing4ti2.converted sing4ti2."+fileending));
+   }
+
    string ergstr = "intvec erglist = " + s + "0;";
    execute(ergstr);
-
+  
    //   print(erglist);
-
+  
    int Rnc = erglist[1];
    int Rnr = erglist[2];
-
+   
    intmat R[Rnr][Rnc];
 
@@ -1910 +2981 @@
      }
    }
+
+
 
    return (R);
@@ -1931 +3004 @@
       1, 1, 0,  0, -1,  0,-1, 0, 0;
    hilbert4ti2intmat(M);
-   hermite(M);
+   hermiteNormalForm(M);
 }
 
@@ -1962 +3035 @@
 
 /******************************************************/
-proc mDegBasis(intvec d)
+proc mDegBasis(intvec d) 
 "
 USAGE: multidegree d
@@ -1991 +3064 @@
   }
 
-  intmat T = getTorsion(R);
+  intmat T = getLattice(R);
 
   if( isFreeRepresented() )
@@ -2006 +3079 @@
 
     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];
+    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 ) ); //
+    intmat K = leftKernelZ(( AA ) ); // 
 
     //      print(K);
@@ -2022 +3095 @@
     //      "!";
 
-    intmat B = hilbert4ti2intmat(transpose(KK), 1);
+    intmat B = hilbert4ti2intmat(transpose(KK), 1);  
 
     //      "!";      print(B);
@@ -2033 +3106 @@
 
 
-  int i;
+  int i;  
   int nnr = nrows(B);
   int nnc = ncols(B);
@@ -2088 +3161 @@
 
   intmat g1[2][2]=1,0,0,1;
-  intmat t[2][1]=2,0;
+  intmat l[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);
@@ -2098 +3171 @@
   v3;
   mDegBasis(v3);
-
-  setBaseMultigrading(g1,t);
+  
+  setBaseMultigrading(g1,l);
   mDegBasis(v1);
   setBaseMultigrading(g2);
   mDegBasis(v2);
-
+  
   intmat M[2][2] = 1, -1, -1, 1;
   intvec d = -2, 2;
@@ -2117 +3190 @@
   intmat M[2][3] = 1, -2, 1,     1, 1, 0;
 
-  intmat T[2][1] = 0, 2;
+  intmat L[2][1] = 0, 2;
 
   intvec d = 4, 1;
 
-  setBaseMultigrading(M, T);
+  setBaseMultigrading(M, L);
 
   mDegBasis(d);
@@ -2162 +3235 @@
 
 proc mDegSyzygy(def I)
-"USAGE: mDegSyzygy(I); I is a poly/vector/ideal/module
+"USAGE: mDegSyzygy(I); I is a ideal or a module
 PURPOSE: computes the multigraded syzygy module of I
 RETURNS: module, the syzygy module of I
@@ -2169 +3242 @@
 "
 {
-  if( isHomogenous(I, "checkGens") == 0)
-  {
-    ERROR ("Sorry: inhomogenous input!");
-  }
+  if( isHomogeneous(I, "checkGens") == 0) 
+  { 
+    ERROR ("Sorry: inhomogeneous input!"); 
+  } 
   module S = syz(I);
   S = setModuleGrading(S, mDeg(I));
@@ -2180 +3253 @@
 {
   "EXAMPLE:"; echo=2;
-
+  
 
   ring r = 0,(x,y,z,w),dp;
@@ -2188 +3261 @@
   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);
+  isHomogeneous(M);
   "Multidegrees: "; print(mDeg(M));
   // Let's compute syzygies!
@@ -2203 +3276 @@
   "Multidegrees: "; print(mDeg(S));
 
-  isHomogenous(S);
+  isHomogeneous(S);
+}
+
+
+
+proc mDegModulo(def I, def J)
+"USAGE: mDegModulo(I); I, J are ideals or modules
+PURPOSE: computes the multigraded 'modulo' module of I and J
+RETURNS: module, see 'modulo' command
+NOTE: I and J should have the same multigrading, and their 
+generators must be multigraded homogeneous
+EXAMPLE: example mDegModulo; shows an example
+"
+{
+  if( (isHomogeneous(I, "checkGens") == 0) or (isHomogeneous(J, "checkGens") == 0) )
+  { 
+    ERROR ("Sorry: inhomogeneous input!"); 
+  } 
+  module K = modulo(I, J);
+  K = setModuleGrading(K, mDeg(I));
+  return (K);
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+
+  ring r = 0,(x,y,z),dp;
+  intmat MM[2][3]=
+    -1,1,1,
+     0,1,3;
+  setBaseMultigrading(MM);
+
+  ideal h1 = x, y, z;
+  ideal h2 = x;
+
+  "Multidegrees: "; print(mDeg(h1));
+   
+  // Let's compute modulo(h1, h2):
+  def K = mDegModulo(h1, h2); K;
+
+  "Module Units Multigrading: "; print( getModuleGrading(K) );
+  "Multidegrees: "; print(mDeg(K));
+
+  isHomogeneous(K);
 }
 
@@ -2210 +3326 @@
 "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
+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!");
-  }
+  if( isHomogeneous(I) == 0) 
+  { 
+    ERROR ("Sorry: inhomogeneous input!"); 
+  } 
 
   def S = groebner(I);
-
+  
   if( typeof(I) == "module" or typeof(I) == "vector" )
   {
-    S = setModuleGrading(S, getModuleGrading(I));
+    S = setModuleGrading(S, getModuleGrading(I));      
   }
 
@@ -2241 +3357 @@
   setBaseMultigrading(MM);
 
-
+  
   module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
-
-
+  
+ 
   intmat v[2][nrows(M)]=
     1,
     0;
-
+  
   M = setModuleGrading(M, v);
 
@@ -2258 +3374 @@
   "Multidegrees: "; print(mDeg(M));
 
-  isHomogenous(M);
+  isHomogeneous(M);
 
   /////////////////////////////////////////////////////////////////////////////
@@ -2266 +3382 @@
   "Multidegrees: "; print(mDeg(S));
 
-  isHomogenous(S);
+  isHomogeneous(S);
 
   /////////////////////////////////////////////////////////////////////////////
@@ -2274 +3390 @@
   "Multidegrees: "; print(mDeg(S));
 
-  isHomogenous(S);
+  isHomogeneous(S);
 }
 
@@ -2281 +3397 @@
 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
+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.
+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!");
-  }
+  if( isHomogeneous(I, "checkGens") == 0) 
+  { 
+    ERROR ("Sorry: inhomogeneous input!"); 
+  } 
 
   def R = res(I, ll, #); list L = R; int l = size(L);
-
+  def V = getModuleGrading(I);
   if( (typeof(I) == "module") or (typeof(I) == "vector") )
   {
-    L[1] = setModuleGrading(L[1], getModuleGrading(I));
-  }
-
-  int i;
+    L[1] = setModuleGrading(L[1], V);
+  }
+
+  int i; 
   for( i = 2; i <= l; i++ )
   {
@@ -2313 +3429 @@
     }
   }
-
+ 
   return (L);
 
-
+  
 }
 example
 {
   "EXAMPLE:"; echo=2;
-
+  
   ring r = 0,(x,y,z,w),dp;
 
@@ -2330 +3446 @@
   setBaseMultigrading(M);
 
-
+  
   module m= ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
-
-  isHomogenous(ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3), "checkGens");
-
+  
+  isHomogeneous(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);
 
@@ -2374 +3490 @@
 
   /////////////////////////////////////////////////////////////////////////////
-
+  
   L = mDegResolution(maxideal(1), 0, 1);
 
@@ -2384 +3500 @@
     "Multigrading: "; print(mDeg(L[j]));
   }
-
+  
   kill v;
-
+  
 
   def h = hilbertSeries(m);
@@ -2393 +3509 @@
   numerator1;
   factorize(numerator1);
-
+  
   denominator1;
   factorize(denominator1);
@@ -2407 +3523 @@
 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.
+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
+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)");
-  }
-
+    "Things might happen, since we are not  free.";
+    //ERROR("SORRY: ONLY TORSION-FREE CASE (POSITIVE GRADING)");
+  }
+   
   int i, j, k, v;
 
   intmat M = getVariableWeights();
-
+  
   int cc = ncols(M);
   int n = nrows(M);
@@ -2437 +3554 @@
 
   int l = size(RES);
-
+  
   list L; L[l + 1] = 0;
 
@@ -2444 +3561 @@
     intmat zeros[n][1];
     L[1] = zeros;
-  }
+  } 
   else
   {
@@ -2454 +3571 @@
     L[j + 1] = mDeg(RES[j]);
   }
-
+  
   l++;
 
   ring R = 0,(t_(1..n),s_(1..n)),dp;
-
-  ideal units;
+  
+  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
@@ -2475 +3592 @@
   poly monom, summand, numerator;
   poly denominator = 1;
-
+  
   for( i = 1; i <= cc; i++)
   {
     monom = 1;
-
     for( k = 1; k <= n; k++)
     {
@@ -2487 +3603 @@
       {
         monom = monom * (var(k)^(v));
-      }
+      } 
       else
       {
@@ -2493 +3609 @@
       }
     }
-
+    
     if( monom == 1)
     {
@@ -2501 +3617 @@
     denominator = denominator * (1 - monom);
   }
-
-  for( j = 1; j<= l; j++)
+  
+  for( j = 1; j<= l; j++) 
   {
     summand = 0;
@@ -2516 +3632 @@
         {
           monom = monom * (var(k)^v);
-        }
+        } 
         else
         {
@@ -2526 +3642 @@
     numerator = numerator - (-1)^j * summand;
   }
-
+  
   if( denominator == 0 )
   {
     ERROR("Multigrading not positive.");
-  }
-
+  } 
+  
   poly denominator1 = denominator;
   poly numerator1 = numerator;
@@ -2565 +3681 @@
   "The s_(i)-variables are defined to be the inverse of the t_(i)-variables.";
   " ------------ ";
-
+  
   return(Q);
 }
@@ -2571 +3687 @@
 {
   "EXAMPLE:"; echo=2;
-
+  
   ring r = 0,(x,y,z,w),dp;
   intmat g[2][4]=
@@ -2577 +3693 @@
     0,1,3,4;
   setBaseMultigrading(g);
-
+  
   module M = ideal(xw-yz, x2z-y3, xz2-y2w, yw2-z3);
   intmat V[2][1]=
@@ -2589 +3705 @@
   factorize(numerator2);
   factorize(denominator2);
-
+  
   kill g, h; setring r;
 
@@ -2595 +3711 @@
     1,2,3,4,
     0,0,5,8;
-
+  
   setBaseMultigrading(g);
-
+  
   ideal I = x^2, y, z^3;
   I = std(I);
@@ -2612 +3728 @@
   mDeg(I);
   def h = hilbertSeries(I); setring h;
-
+  
   factorize(numerator2);
   factorize(denominator2);
@@ -2619 +3735 @@
   ////////////////////////////////////////////////
   ring R = 0,(x,y,z),dp;
-  intmat W[2][3] =
+  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);
@@ -2633 +3749 @@
   ////////////////////////////////////////////////
   ring R = 0,(x,y,z,a,b,c),dp;
-  intmat W[2][6] =
+  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] =
+  intmat W[1][4] = 
      1,1, 1,1;
   setBaseMultigrading(W);
@@ -2655 +3771 @@
 
   hilb(std(I));
-
+  
   def h = hilbertSeries(I); setring h;
-
+  
   numerator1;
   denominator1;
@@ -2663 +3779 @@
   factorize(numerator2);
   factorize(denominator2);
-
+  
 
   kill h;
@@ -2672 +3788 @@
 
   hilb(std(I2));
-
+  
   def h = hilbertSeries(I2); setring h;
 
@@ -2682 +3798 @@
   ////////////////////////////////////////////////
   setring R;
-
+  
   W = 2,2,2,2;
-
+  
   setBaseMultigrading(W);
 
@@ -2694 +3810 @@
 
   kill w;
-
+  
 
   def h = hilbertSeries(I2); setring h;
 
-
+  
   numerator1; denominator1;
   kill h;
 
-
+  
   kill R, W;
 
@@ -2716 +3832 @@
 
   hilb(std(I));
-
+  
   def h = hilbertSeries(I); setring h;
 
@@ -2732 +3848 @@
 
   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;
@@ -2737 +3879 @@
   setring R;
 
-  I = x^5; I;
-
-  hilb(std(I));
-  hilb(std(I), 1);
-
-  def h = hilbertSeries(I); setring h;
+  module M = 1;
+
+  M = setModuleGrading(M, W);
+
+  
+  hilb(std(M));
+
+  def h = hilbertSeries(M); setring h;
 
   numerator1; denominator1;
-
 
   kill h;
@@ -2751 +3894 @@
   setring R;
 
-  I = x^10; I;
-
-  hilb(std(I));
-
-  def h = hilbertSeries(I); setring h;
+  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;
@@ -2761 +3908 @@
   kill h;
   ////////////////////////////////////////////////
-  setring R;
-
-  module M = 1;
-
-  M = setModuleGrading(M, W);
-
-
-  hilb(std(M));
-
-  def h = hilbertSeries(M); setring 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;
@@ -2776 +3926 @@
   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;
-
-
+  setring R; 
+
+  
   module S = M + N;
-
+  
   S = setModuleGrading(S, V);
 
@@ -2830 +3946 @@
 }
 
-
+proc evalHilbertSeries(def h, intvec v)
+"
+evaluate hilbert series h by substibuting v[i] for t_(i) (1/v[i] for s_(i))
+return: int (h(v))
+"
+{
+  if( 2*size(v) != nvars(h) ) 
+  {
+    ERROR("Wrong input/size!");
+  }
+
+  setring h;
+
+  if( defined(numerator2) and defined(denominator2) )
+  {
+    poly n = numerator2; poly d = denominator2;
+  } else
+  {
+    poly n = numerator1; poly d = denominator1;
+  }
+
+  int N = size(v);
+  int i; number k;
+  ideal V;
+
+  for( i = N; i > 0; i -- )
+  {
+    k = v[i];
+    V[i] = var(i) - k;
+  }
+   
+  V = groebner(V);
+   
+  n = NF(n, V); 
+  d = NF(d, V); 
+
+  n;
+  d;
+
+  if( d == 0 )
+  {
+    ERROR("Sorry: denominator is zero!");
+  }
+  
+  if( n == 0 )
+  {
+    return (0);
+  }
+
+  poly g = gcd(n, d);
+  
+  if( g != leadcoef(g) )
+  {
+    n = n / g;
+    d = d / g;
+  }
+
+  n;
+  d;
+   
+   
+  for( i = N; i > 0; i -- )
+  {
+    "i: ", i;
+    n;
+    d;
+     
+    k = v[i];
+    k;
+     
+    n = subst(n, var(i), k);
+    d = subst(d, var(i), k);
+    
+    if( k != 0 )
+    {
+      k = 1/k;
+      n = subst(n, var(N+i), k);
+      d = subst(d, var(N+i), k);
+    }
+  }
+
+  n;
+  d;
+
+  if( d == 0 )
+  {
+    ERROR("Sorry: denominator is zero!");
+  }
+  
+  if( n == 0 )
+  {
+    return (0);
+  }
+
+  poly g = gcd(n, d);
+  
+  if( g != leadcoef(g) )
+  {
+    n = n / g;
+    d = d / g;
+  }
+
+  n;
+  d;
+   
+  if( n != leadcoef(n) || d != leadcoef(d) )
+  {
+    ERROR("Sorry cannot completely evaluate. Partial result: (" + string(n) + ")/(" + string(d) + ")");
+  }
+
+  n;
+  d;
+
+  return (leadcoef(n)/leadcoef(d));
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+
+  // TODO!
+
+}
+
+
+proc isPositive()
+"USAGE: isPositive()
+PURPOSE: Computes whether the multigrading of the ring is positive.
+For computation theorem 8.6 of the Miller/Sturmfels book is used.
+RETURNS: true if the multigrading is positive
+EXAMPLE: example isPositive; shows an example
+"
+{
+ideal I = mDegBasis(0);
+ideal J = attrib(I,"ZeroPart");
+/* 
+I am not quite sure what this ZeroPart is anymore. I thought it
+should contain all monomials of degree 0, but then apparently 1 should
+be contained. It makes sense to exclude 1, but was this also the intention?
+*/
+return(J==0);
+}
+example
+{
+  echo = 2; printlevel = 3;
+  ring r = 0,(x,y),dp;
+  intmat A[1][2]=-1,1;
+  setBaseMultigrading(A);
+  isPositive();
+  
+  intmat A[1][2]=1,1;
+  setBaseMultigrading(A);
+  isPositive(A);
+}
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -2836 +4104 @@
 proc testMultigradingLib ()
 {
-  echo = 2; printlevel = 3;
   example setBaseMultigrading;
   example setModuleGrading;
 
   example getVariableWeights;
-  example getTorsion;
+  example getLattice;
+  example getGradingGroup;
   example getModuleGrading;
 
@@ -2849 +4117 @@
 
 
-  example hermite;
-  example isHomogenous;
+  example hermiteNormalForm;
+  example isHomogeneous;
   example isTorsionFree;
   example pushForward;
-  example defineHomogenous;
+  example defineHomogeneous;
 
   example equalMDeg;
-  example isTorsionElement;
+  example isZeroElement;
 
   example mDegResolution;
-
-  "// ******************* example hilbertSeries ************************//";
+  
+  "// ******************* example hilbertSeries ************************//"; 
   example hilbertSeries;
 
@@ -2867 +4135 @@
 
   "The End!";
-
-}
+}
+
+
+proc mDegTruncate(def M, intvec md)
+{
+  "d: ";
+  print(md);
+  
+  "M: ";
+  module LL = M; // + L for d+1
+  LL;
+  print(mDeg(LL));
+
+
+  intmat V = getModuleGrading(M); 
+  intvec vi;
+  int s = nrows(M); 
+  int r = nrows(V);
+  int i;
+  module L; def B;
+  for (i=s; i>0; i--) 
+  {
+    "comp: ", i;
+    vi = V[1..r, i];
+    "v[i]: "; vi;
+
+    B = mDegBasis(md - vi); // ZeroPart is always the same...
+    "B: "; B;
+
+    L = L, B*gen(i);
+  }
+  L = simplify(L, 2);
+  L = setModuleGrading(L,V);
+
+  "L: "; L;
+  print(mDeg(L));
+
+  L = mDegModulo(L, LL);
+  L = mDegGroebner(L);
+//  L = minbase(prune(L));
+
+  "??????????";
+  print(L);
+  print(mDeg(L));
+   
+  V = getModuleGrading(L);
+
+  // take out other degrees
+  for(i = ncols(L); i > 0; i-- )
+  {
+    if( !equalMDeg( mDeg(getGradedGenerator(L, i)), md ) )
+    {
+      L[i] = 0;
+    }
+  }
+  
+  L = simplify(L, 2);
+  L = setModuleGrading(L, V);
+  print(L);
+  print(mDeg(L));
+  
+  return(L);
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+
+  // TODO!
+  ring r = 32003, (x,y), dp;
+  
+  intmat M[2][2] = 
+    1, 0, 
+    0, 1;
+
+  setBaseMultigrading(M);
+
+  intmat V[2][1] = 
+    0, 
+    0;
+  
+  "X:";
+  module h1 = x;
+  h1 = setModuleGrading(h1, V);
+  mDegTruncate(h1, mDeg(x));
+  mDegTruncate(h1, mDeg(y));
+
+  "XY:";
+  module h2 = ideal(x, y);
+  h2 = setModuleGrading(h2, V);
+  mDegTruncate(h2, mDeg(x));
+  mDegTruncate(h2, mDeg(y));
+  mDegTruncate(h2, mDeg(xy));
+}
+
+
+/******************************************************/
+/* Some functions on lattices. 
+TODO Tuebingen: - add functionality (see wiki) and 
+- adjust them to work for groups as well.*/
+/******************************************************/
+
+
+
+/******************************************************/
+proc imageLattice(intmat Q, intmat L)
+"USAGE: imageLattice(Q,L); Q and L are of type intmat
+PURPOSE: compute an integral basis for the image of the 
+lattice L under the homomorphism of lattices Q.
+RETURN: intmat
+EXAMPLE: example imageLattice; shows an example
+"
+{
+  intmat Mul = Q*L;
+  intmat L = latticeBasis(Mul);
+
+  return(L);
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+   
+  intmat Q[2][3] =
+    1,2,3,
+    3,2,1;
+
+  intmat L[3][2] =
+    1,4,
+    2,5,
+    3,6;
+
+  // should be a 2x2 matrix with columns
+  // [2,-14], [0,36]
+  imageLattice(Q,L);
+
+}
+
+/******************************************************/
+proc intRank(intmat A)
+"
+USAGE: intRank(A); intmat A
+PURPOSE: compute the rank of the integral matrix A 
+by computing a hermite normalform.
+RETURNS: int
+EXAMPLE: example intRank; shows an example
+"
+{
+  intmat B = hermiteNormalForm(A);
+
+  // get number of zero columns
+  int nzerocols = 0;
+  int j;
+  int i;
+  int iszero;
+  for ( j = 1; j <= ncols(B); j++ )
+  {
+    iszero = 1;
+    
+    for ( i = 1; i <= nrows(B); i++ )
+    {
+      if ( B[i,j] != 0 ) 
+      {
+        iszero = 0;
+        break;
+      }
+    }
+    
+    if ( iszero == 1 )
+    {
+      nzerocols++;
+    }
+  }
+
+  // get number of zero rows
+  int nzerorows = 0;
+
+  for ( i = 1; i <= nrows(B); i++ )
+  {
+    iszero = 1;
+    
+    for ( j = 1; j <= ncols(B); j++ )
+    {
+      if ( B[i,j] != 0 ) 
+      {
+        iszero = 0;
+        break;
+      }
+    }
+    
+    if ( iszero == 1 )
+    {
+      nzerorows++;
+    }
+  }
+
+  int r = nrows(B) - nzerorows;
+
+  if ( ncols(B) - nzerocols < r ) 
+  {
+    r = ncols(B) - nzerocols;
+  }
+  
+  return(r);
+}
+example
+{
+
+  intmat A[3][4] =
+    1,0,1,0,
+    1,2,0,0,
+    0,0,0,0;
+
+  int r = intRank(A);
+
+  print(A);
+  print(r); // Should be 2
+
+  kill A;
+
+}
+
+/*****************************************************/
+
+proc isSublattice(intmat L, intmat S)
+"USAGE: isSublattice(L, S); L, S are of tpye intmat
+PURPOSE: checks whether the lattice created by L is a 
+sublattice of the lattice created by S.
+The procedure checks whether each generator of L is 
+contained in S.
+RETURN: 0 if false, 1 if true
+EXAMPLE: example isSublattice; shows an example
+"
+{
+  int a,b,g,i,j,k;
+  intmat Ker;
+    
+  // check whether each column v of L is contained in
+  // the lattice generated by S
+  for ( i = 1; i <= ncols(L); i++ )
+  {
+    
+    // v is the i-th column of L
+    intvec v;
+     for ( j = 1; j <= nrows(L); j++ )
+    {
+      v[j] = L[j,i];
+    }
+
+    // concatenate B = [S,v]
+    intmat B[nrows(L)][ncols(S) + 1];
+
+    for ( a = 1; a <= nrows(S); a++ )
+    {
+      for ( b = 1; b <= ncols(S); b++ )
+      {
+        B[a,b] = S[a,b];
+      }
+    }
+
+    for ( a = 1; a <= size(v); a++ )
+    {
+      B[a,ncols(B)] = v[a];
+    }
+
+    
+    // check gcd
+    Ker = kernelLattice(B);
+    k = nrows(Ker);
+    list R; // R is the last row
+
+    for ( j = 1; j <= ncols(Ker); j++ )
+    {
+      R[j] = Ker[k,j];
+    }
+
+    g = R[1];
+    
+    for ( j = 2; j <= size(R); j++ )
+    {
+      g = gcd(g,R[j]);
+    }
+
+    if ( g != 1 ) 
+    {
+      return(0);
+    }
+
+    kill B, v, R;
+
+  }
+
+  return(1);
+}
+example
+{
+  "EXAMPLE:"; echo=2;
+   
+  //ring R = 0,(x,y),dp;
+  intmat S2[3][2]=
+    0, 2, 3,
+    0, 1, 1,
+    3, 0, 2;
+
+  intmat S1[3][3]=
+    0, 6,
+    0, 2,
+    3, 4;
+
+  isSublattice(S1,S2); // Yes!
+
+  intmat S3[3][1] =
+    0,
+    0,
+    1;
+
+  not(isSublattice(S3,S2)); // Yes!
+
+}
+
+/******************************************************/
+
+proc latticeBasis(intmat B)
+"USAGE: latticeBasis(B); intmat B
+PURPOSE: compute an integral basis for the lattice defined by 
+the columns of B.
+RETURNS: intmat
+EXAMPLE: example latticeBasis; shows an example
+"
+{
+  int n = ncols(B);
+  int r = intRank(B); 
+ 
+  if ( r == 0 ) 
+  {
+    intmat H[nrows(B)][1];
+    int j;
+
+    for ( j = 1; j <= nrows(B); j++ )
+    {
+      H[j,1] = 0;   
+    }
+  }
+  else
+  {
+    intmat H = hermiteNormalForm(B);;
+
+    if (r < n) 
+    {
+      // delete columns r+1 to n
+      // should be identical with the function
+      // H = submat(H,1..nrows(H),1..r); 
+      // for matrices
+      intmat Hdel[nrows(H)][r];
+      int k;
+      int m;
+      
+      for ( k = 1; k <= nrows(H); k++ )
+      {
+        for ( m = 1; m <= r; m++ )
+        {
+          Hdel[k,m] = H[k,m];
+        }
+      }
+
+      H = Hdel;      
+    }
+  }
+  
+  return(H);  
+} 
+example
+{
+  "EXAMPLE:"; echo=2;
+  
+  intmat L[3][3] =
+    1,4,8,
+    2,5,10,
+    3,6,12;
+
+  intmat B = latticeBasis(B); // should be a matrix with columns [1,2,3] and [0,3,6]
+
+  kill B,L;
+}
+
+/******************************************************/
+
+proc kernelLattice(def P)
+"
+USAGE: kernelLattice(P); intmat P
+PURPOSE: compute a integral basis for the kernel of the
+homomorphism of lattices defined by the intmat P.
+RETURNS: intmat
+EXAMPLE: example kernelLattice; shows an example
+"
+{
+  int n = ncols(P);
+  int r = intRank(P);
+
+  if ( r == 0 )
+  {
+    intmat U = unitMatrix(n);
+  }
+  else
+  {
+    if ( r == n ) 
+    {
+      intmat U[n][1];  // now all entries are zero.
+    }
+    else
+    {
+      list L = hermiteNormalForm(P, "transform"); //hermite(P, "transform");  // now, Hermite = L[1] = A*L[2]
+      intmat U = L[2];
+
+      // delete columns 1 to r
+      // should be identical with the function
+      // U = submat(U,1..nrows(U),r+1..); 
+      // for matrices
+      intmat Udel[nrows(U)][ncols(U) - r];
+      int k;
+      int m;
+      
+      for ( k = 1; k <= nrows(U); k++ )
+      {
+        for ( m = r + 1; m <= ncols(U); m++ )
+        {
+          Udel[k,m - r] = U[k,m];
+        }
+      }
+
+      U = Udel;  
+
+    }
+  }
+
+  return(U);
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+
+  intmat LL[3][4] = 
+    1,4,7,10,
+    2,5,8,11,
+    3,6,9,12;
+
+  // should be a 4x2 matrix with colums
+  // [-1,2,-1,0],[2,-3,0,1]
+  intmat B = kernelLattice(LL);
+
+  print(B);
+
+  kill B;
+
+}
+
+/*****************************************************/
+
+proc preimageLattice(def P, def B)
+"
+USAGE: preimageLattice(P, B); intmat P, intmat B
+PURPOSE: compute an integral basis for the preimage of B under
+ the homomorphism of lattices defined by the intmat P.
+RETURNS: intmat
+EXAMPLE: example preimageLattice; shows an example
+"
+{
+  // concatenate matrices: Con = [P,-B]
+  intmat Con[nrows(P)][ncols(P) + ncols(B)];
+  int i;
+  int j;
+
+  for ( i = 1; i <= nrows(Con); i++ )
+  {
+    for ( j = 1; j <= ncols(P); j++ ) // P first
+    {
+      Con[i,j] = P[i,j];
+    }
+  }
+
+  for ( i = 1; i <= nrows(Con); i++ )
+  {
+    for ( j = 1; j <= ncols(B); j++ ) // now -B
+    {
+      Con[i,ncols(P) + j] = - B[i,j];
+    }
+  }
+
+
+  print(Con);
+
+  intmat L = kernelLattice(Con);
+
+  print(L);
+  print(ncols(P));
+  print(ncols(L));
+
+  // delete rows ncols(P)+1 to nrows(L) out of L
+  intmat Del[ncols(P)][ncols(L)];
+  int k;
+  int m;
+      
+  for ( k = 1; k <= nrows(Del); k++ )
+  {
+    for ( m = 1; m <= ncols(Del); m++ )
+    {
+      Del[k,m] = L[k,m];
+    }
+  }
+  
+  L = latticeBasis(Del);
+
+  return(L); 
+
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+
+  intmat P[2][3] = 
+    2,6,10,
+    4,8,12;
+
+  intmat B[2][1] =
+    1,
+    0;
+
+  // should be a 2x2 matrix with columns
+  // [1,1,-1] and [-2,1,0]
+  intmat L = preimageLattice(P,B);
+
+  kill B, P, L;
+
+}
+
+/******************************************************/
+proc isPrimitiveSublattice(intmat A);
+"USAGE: isPrimitiveSublattice(A); intmat A
+PURPOSE: check whether the given set of integral vectors in ZZ^m, 
+i.e. the columns of A, generate a primitive sublattice in ZZ^m 
+(a direct summand of ZZ^m). 
+RETURNS: int, where 0 is false and 1 is true.
+EXAMPLE: example isPrimitiveSublattice; shows an example
+"
+{
+  intmat B = smithNormalForm(A);
+  int r = intRank(B);
+  
+  if ( r == 0 ) 
+  {
+    return(1);
+  }
+
+  if ( 1 < B[r,r] )
+  {
+    return(0);
+  }
+
+  return(1);
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+  
+  intmat A[3][2] =
+    1,4,
+    2,5,
+    3,6;
+
+  // should be 0
+  int b = isPrimitiveSublattice(A);
+ 
+  kill A,b;
+}
+
+/******************************************************/
+proc isIntegralSurjective(intmat P);
+"USAGE: isIntegralSurjective(P); intmat P
+PURPOSE: test whether the given linear map P of lattices is 
+surjective.
+RETURNS: int, where 0 is false and 1 is true.
+EXAMPLE: example isIntegralSurjective; shows an example
+"
+{
+  int r = intRank(P);
+  
+  if ( r < nrows(P) )
+  {
+    return(0);
+  }
+
+  if ( isPrimitiveSublattice(A) == 1 )
+  {
+    return(1);
+  }
+
+  return(0);
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+  
+  intmat A[3][2] =
+    1,3,5,
+    2,4,6;
+    
+  // should be 0
+  int b = isIntegralSurjective(A);
+  print(b);
+ 
+  kill A,b;
+}
+
+/******************************************************/
+proc projectLattice(intmat B)
+"USAGE: projectLattice(B); intmat B
+PURPOSE: A set of vectors in ZZ^m is given as the columns of B. 
+Then this function provides a linear map ZZ^m --> ZZ^n 
+having the primitive span of B its kernel.
+RETURNS: intmat
+EXAMPLE: example projectLattice; shows an example
+"
+{
+  int n = nrows(B);
+  int r = intRank(B);
+
+  if ( r == 0 )
+  {
+    intmat U = unitMatrix(n);
+  }
+  else
+  {
+    if ( r == n ) 
+    {
+      intmat U[n][1]; // U now is the n-dim zero-vector
+    }
+    else
+    {
+      // we want a matrix with column operations so we transpose
+      list L = hermiteNormalForm(B, "transform"); //hermite(transpose(B), "transform"); 
+      intmat U = transpose(L[2]); 
+
+      // delete rows 1 to r
+      intmat Udel[nrows(U) - r][ncols(U)];
+      int k;
+      int m;
+      
+      for ( k = 1; k <= nrows(U) - r ; k++ )
+      {
+        for ( m = 1; m <= ncols(U); m++ )
+        {
+          Udel[k,m] = U[k + r,m];
+        }
+      }
+
+      U = Udel;  
+      
+    }
+  }
+  
+  return(U);
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+  
+  intmat B[4][2] = 
+    1,5,
+    2,6,
+    3,7,
+    4,8;
+  
+  // should result in a (2x4)-matrix with columns
+  // [-1, 2], [2, −3], [-1, 0] and [0, 1].
+  intmat U = projectLattice(B);
+  
+}
+
+/******************************************************/
+proc intersectLattices(intmat A, intmat B)
+"USAGE: intersectLattices(A, B); intmat A, intmat B
+PURPOSE: compute an integral basis for the intersection of the 
+lattices A and B.
+RETURNS: intmat
+EXAMPLE: example intersectLattices; shows an example
+"
+{
+  // concatenate matrices: Con = [A,-B]
+  intmat Con[nrows(A)][ncols(A) + ncols(B)];
+  int i;
+  int j;
+
+  for ( i = 1; i <= nrows(Con); i++ )
+  {
+    for ( j = 1; j <= ncols(A); j++ ) // A first
+    {
+      Con[i,j] = A[i,j];
+    }
+  }
+
+  for ( i = 1; i <= nrows(Con); i++ )
+  {
+    for ( j = 1; j <= ncols(B); j++ ) // now -B
+    {
+      Con[i,ncols(A) + j] = - B[i,j];
+    }
+  }
+
+  intmat K = kernelLattice(Con);
+
+  // delete all rows in K from ncols(A)+1 onwards
+  intmat Bas[ncols(A)][ncols(K)];
+  
+  for ( i = 1; i <= nrows(Bas); i++ )
+  {
+    for ( j = 1; j <= ncols(Bas); j++ ) 
+    {
+      Bas[i,j] = K[i,j];
+    }
+  }
+
+  // take product in order to obtain the intersection
+  intmat S = A * Bas;
+  intmat Cut = hermiteNormalForm(S); //hermite(S);
+  int r = intRank(Cut);
+
+  if ( r == 0 ) 
+  {
+    intmat Cutdel[nrows(Cut)][1]; // is now the zero-vector
+
+    Cut = Cutdel;
+  }
+  else
+  {
+    // delte columns from r+1 onwards
+    intmat Cutdel[nrows(Cut)][r];
+
+    for ( i = 1; i <= nrows(Cutdel); i++ )
+    {
+      for ( j = 1; j <= r; j++ ) 
+      {
+        Cutdel[i,j] = Cut[i,j];
+      }
+    }
+
+    Cut = Cutdel;
+  }
+
+  return(Cut);
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+  
+  intmat A[3][2] = 
+    1,4,
+    2,5,
+    3,6;
+
+  intmat B[3][2] = 
+    6,9,
+    7,10,
+    8,11;
+  
+  // should result in a (2x4)-matrix with columns
+  // [3, 3, 3], [0, 3, 6]
+  intmat U = intersectLattices(A,B);
+  
+}
+
+proc intInverse(intmat A);
+"USAGE: intInverse(A); intmat A
+PURPOSE: compute the integral inverse of the intmat A. 
+If det(A) is neither 1 nor -1 an error is returned.
+RETURNS: intmat
+EXAMPLE: example intInverse; shows an example
+"
+{
+  int d = det(A);
+  
+  if ( d * d != 1 ) // is d = 1 or -1? Else: error
+  {
+    ERROR("determinant of the given intmat has to be 1 or -1.");
+  }
+  
+  int c;
+  int i,j;
+  intmat C[nrows(A)][ncols(A)];
+  intmat Ad;
+  int s;
+
+  for ( i = 1; i <= nrows(C); i++ )
+  {
+    for ( j = 1; j <= ncols(C); j++ )
+    {
+      Ad = intAdjoint(A,i,j);
+      s = 1;
+    
+      if ( ((i + j) % 2) > 0 )
+      {
+        s = -1;
+      }
+
+      C[i,j] = d * s * intDet(Ad); // mult by d is equal to div by det
+    }
+  }
+
+  C = transpose(C);
+
+  return(C);
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+
+  intmat A[3][3] =
+    1,1,3,
+    3,2,0,
+    0,0,1;
+
+  intmat B = intInverse(A);
+
+  // should be the unit matrix
+  print(A * B);
+
+  kill A,B;
+}
+
+
+/******************************************************/
+proc intAdjoint(intmat A, int indrow, int indcol)
+"USAGE: intAdjoint(A); intmat A
+PURPOSE: return the matrix where the given row and column are deleted.
+RETURNS: intmat
+EXAMPLE: example intAdjoint; shows an example
+"
+{
+  int n = nrows(A);
+  int m = ncols(A);
+  int i, j;
+  intmat B[n - 1][m - 1];
+  int a, b;
+
+  for ( i = 1; i < indrow; i++ )
+  {
+    for ( j = 1; j < indcol; j++ )
+    {
+      B[i,j] = A[i,j];
+    }
+    for ( j = indcol + 1; j <= ncols(A); j++ )
+    {
+      B[i,j - 1] = A[i,j];
+    }
+  }
+
+  for ( i = indrow + 1; i <= nrows(A); i++ )
+  {
+    for ( j = 1; j < indcol; j++ )
+    {
+      B[i - 1,j] = A[i,j];
+    }
+    for ( j = indcol+1; j <= ncols(A); j++ )
+    {
+      B[i - 1,j - 1] = A[i,j];
+    }
+  }
+
+  return(B);
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+  
+  intmat A[2][3] =
+    1,3,5,
+    2,4,6;
+
+  intmat B = intAdjoint(A,2,2);
+  print(B);
+
+  kill A,B;
+}
+
+/******************************************************/
+proc integralSection(intmat P);
+"USAGE: integralSection(P); intmat P
+PURPOSE: for a given linear surjective map P of lattices
+ this procedure returns an integral section of P.
+RETURNS: intmat
+EXAMPLE: example integralSection; shows an example
+"
+{
+  int m = nrows(P);
+  int n = ncols(P);
+  
+  if ( m == n )
+  {
+    intmat U = intInverse(P); 
+  }
+  else
+  {
+    intmat U = (hermiteNormalForm(P, "transform"))[2];
+    
+    // delete columns m+1 to n
+    intmat Udel[nrows(U)][ncols(U) - (n - m)];
+    int k;
+    int z;
+      
+    for ( k = 1; k <= nrows(U); k++ )
+    {
+      for ( z = 1; z <= m; z++ )
+      {
+        Udel[k,z] = U[k,z];
+      }
+    }
+    
+    U = Udel;  
+  }
+
+  return(U);
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+
+  intmat P[2][4] =
+    1,3,4,6,
+    2,4,5,7;
+
+  // should be a matrix with two columns 
+  // for example: [−2, 1, 0, 0], [3, −3, 0, 1]
+  intmat U = integralSection(P);
+
+  print(U);
+  print(P * U);
+
+  kill U;  
+}
+
+
+
+/******************************************************/
+proc factorgroup(G,H)
+"USAGE: factorgroup(G,H); list G, list H
+PURPOSE: returns a representation of the factor group G mod H using the first isomorphism thm
+RETURNS: list
+EXAMPLE: example factorgroup(G,H); shows an example
+"
+{
+  intmat S1 = G[1];
+  intmat L1 = G[2];
+  intmat S2 = H[1];
+  intmat L2 = H[2];
+
+  // check whether G,H are subgroups of a common group, i.e. whether L1 and L2 span the same lattice  
+  if ( !isSublattice(L1,L2) || !isSublattice(L2,L1))
+  {
+    ERROR("G and H are not subgroups of a common group.");
+  }
+
+  // check whether H is a subgroup of G, i.e. whether S2 is a sublattice of S1+L1
+  intmat B = concatintmat(S1,L1); // check whether this gives the concatinated matrix
+  if ( !isSublattice(S2,B) )
+  {
+    ERROR("H is not a subgroup of G");
+  }
+  // use first isomorphism thm to get the factor group
+  intmat L = concatintmat(L1,S2); // check whether this gives the concatinated matrix
+  list GmodH;
+  GmodH[1]=S1;
+  GmodH[2]=L;
+  return(GmodH);
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+  
+  intmat S1[2][2] =
+    1,0,
+    0,1;
+  intmat L1[2][1] =
+    2,
+    0;
+
+  intmat S2[2][1] = 
+    1,
+    0;
+  intmat L2[2][1] =
+    2,
+    0;
+
+  list G = createGroup(S1,L1);
+  list H = createGroup(S2,L2);
+
+  list N = factorgroup(G,H);
+  print(N);
+
+  kill G,H,N,S1,L1,S2,L2;
+    
+}
+
+/******************************************************/
+proc productgroup(G,H)
+"USAGE: productgroup(G,H); list G, list H
+PURPOSE: returns a representation of the G x H
+RETURNS: list
+EXAMPLE: example productgroup(G,H); shows an example
+"
+{
+  intmat S1 = G[1];
+  intmat L1 = G[2];
+  intmat S2 = H[1];
+  intmat L2 = H[2];
+  intmat OS1[nrows(S1)][ncols(S2)];
+  intmat OS2[nrows(S2)][ncols(S1)];
+  intmat OL1[nrows(L1)][ncols(L2)];
+  intmat OL2[nrows(L2)][ncols(L1)];
+
+  // concatinate matrices to get S
+  intmat A = concatintmat(S1,OS1); 
+  intmat B = concatintmat(OS2,S2); 
+  intmat At = transpose(A);
+  intmat Bt = transpose(B);
+  intmat St = concatintmat(At,Bt);
+  intmat S = transpose(St);
+
+  // concatinate matrices to get L
+  intmat C = concatintmat(L1,OL1); 
+  intmat D = concatintmat(OL2,L2); 
+  intmat Ct = transpose(C);
+  intmat Dt = transpose(D);
+  intmat Lt = concatintmat(Ct,Dt);
+  intmat L = transpose(Lt);
+
+  list GxH;
+  GxH[1]=S;
+  GxH[2]=L;
+  return(GxH);
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+  
+  intmat S1[2][2] =
+    1,0,
+    0,1;
+  intmat L1[2][1] =
+    2,
+    0;
+
+  intmat S2[2][2] = 
+    1,0,
+    0,2;
+  intmat L2[2][1] =
+    0,
+    3;
+
+  list G = createGroup(S1,L1);
+  list H = createGroup(S2,L2);
+
+  list N = productgroup(G,H);
+  print(N);
+
+  kill G,H,N,S1,L1,S2,L2;
+    
+}
+
+/******************************************************/
+proc primitiveSpan(intmat V);
+"USAGE: isIntegralSurjective(V); intmat V
+PURPOSE: compute an integral basis for the minimal primitive
+sublattice that contains the given vectors, i.e. the columns of V.
+RETURNS: int, where 0 is false and 1 is true.
+EXAMPLE: example isIntegralSurjective; shows an example
+"
+{
+  int n = ncols(V);
+  int m = nrows(V);
+  int r = intRank(V);
+
+ 
+  if ( r == 0 )
+  {
+    intmat P[m][1]; //  this is the m-zero-vector now
+  }
+  else
+  {
+    list L = smithNormalForm(V, "transform"); // L = [A,S,B] where S is the smith-NF and S = A*S*B
+    intmat P = intInverse(L[1]); 
+
+    print(L);
+    
+    if ( r < m ) 
+    {
+      // delete columns r+1 to m in P:
+      intmat Pdel[nrows(P)][r];
+      int i,j;
+
+      for ( i = 1; i <= nrows(Pdel); i++ )
+      {
+        for ( j = 1; j <= ncols(Pdel); j++ )
+        {
+          Pdel[i,j] = P[i,j];
+        }
+      }
+
+      P = Pdel;
+    }
+  }
+
+  return(P);
+}
+example
+{
+  "EXAMPLE"; echo = 2;
+  
+  intmat V[2][4] =
+    1,4,
+    2,5,
+    3,6;
+
+  // should return a (4x2)-matrix with columns
+  // [1, 2, 3] and [1, 1, 1] (or similar)
+  intmat R = primitiveSpan(V);
+  print(R);
+
+  kill V,R;
+}
+
+/***********************************************************/