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

// -*- 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$";
category="Combinatorial Commutative Algebra";
info="
LIBRARY:  multigrading.lib          Multigraded Rings

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}
@*        Anna-Lena Winz, anna-lena.winz@math.fu-berlin.de

OVERVIEW: This library allows one to virtually add multigradings to Singular:
grade multivariate polynomial rings with arbitrary (fin. gen. Abelian) groups. 
For more see http://code.google.com/p/convex-singular/wiki/Multigrading
For theoretical references see:
@* E. Miller, B. Sturmfels: 'Combinatorial Commutative Algebra'
and
@* M. Kreuzer, L. Robbiano: 'Computational Commutative Algebra'.

NOTE: 'multiDegBasis' relies on 4ti2 for computing Hilbert Bases.
All groups are finitely generated Abelian

PROCEDURES:
setBaseMultigrading(M,L); attach multiweights/grading group matrices to the basering
getVariableWeights([R]);  get matrix of multidegrees of vars 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

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 P(L)
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);          computes 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.
projectLattice(B);        computes a linear map of lattices having the primitive span of B as its kernel.
intersectLattices(A,B);   computes 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);            computes the integral inverse matrix 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);         computes 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

multiDeg(A);                  compute the multidegree of A
multiDegBasis(d);             compute all monomials of multidegree d
multiDegPartition(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

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

multiDegGroebner(M);          compute the multigraded GB/SB of M
multiDegSyzygy(M);            compute the multigraded syzygies of M
multiDegModulo(I,J);          compute the multigraded 'modulo' module of I and J
multiDegResolution(M,l[,m]);  compute the multigraded resolution of M
multiDegTensor(m,n);          compute the tensor product of multigraded modules m,n
multiDegTor(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
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

lll(A);                   applies LLL(.) of lll.lib which only works for lists on a matrix A

           (parameters in square brackets [] are optional)

KEYWORDS:  multigrading, multidegree, multiweights, multigraded-homogeneous, integral linear algebra
";

/// evalHilbertSeries(h,v);   evaluate hilberts series h by substituting v[i] for t_(i) (too experimentall)

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

LIB "standard.lib"; // for groebner
// LIB "lll.lib"; // for lll_matrix // no need now, right?
LIB "matrix.lib"; // for multiDegTor

/******************************************************/

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(src, Im, 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;

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


/******************************************************/
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 = multiDeg(Im);
//  print(result_degs);

  setring src;

  intmat input_degs = multiDeg(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); // FALSE: there is no such homomorphism!
  // Therefore: the following command should return an error
  // createGradedRingHomomorphism(r, ideal(F), B);

}


proc createQuotientGroup(intmat L)
"USAGE: createGroup(L); L is an integer matrix
PURPOSE: create the group of the form (I+L)/L,
where I is the square identity matrix of size nrows(L) x nrows(L)
NOTE: L specifies relations between free generators of Z^nrows(L)
RETURN: group
EXAMPLE: example createQuotientGroup; shows an example
"
{
  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));
}
example
{
  "EXAMPLE:"; echo=2;

  intmat I[3][3] =
                  1, 0, 0,
                  0, 1, 0,
                  0, 0, 1;

  intmat L[3][2] =
                  1, 1,
                  1, 3,
                  1, 5;


  // The group Z^3 / L can be constructed as follows:

  // shortcut:
  def G = createQuotientGroup(L);
  printGroup(G);

  // the general way:
  def GG = createGroup(I, L); // (I+L)/L
  printGroup(GG);
}

proc createTorsionFreeGroup(intmat S)
"USAGE: createTorsionFreeGroup(S); S is an integer matrix
PURPOSE: create the free subgroup generated by S within the
free Abelian group of rank nrows(S)
RETURN: group
EXAMPLE: example createTorsionFreeGroup; shows an example
"
{
  int r = nrows(S); int i;
  intmat L[r][1] = 0;
  return (createGroup(S,L));
}
example
{
  "EXAMPLE:"; echo=2;

  // ----------- extreme case ------------ //
  intmat S[1][3] =
                  1,  -1, 10;

  // Torsion:
  intmat L[1][1] =
                  0;


  // The free subgroup generated by elements of S within Z^1
  // can be constructed as follows:

  // shortcut:
  def G = createTorsionFreeGroup(S);
  printGroup(G);

  // the general way:
  def GG = createGroup(S, L); // (S+L)/L
  printGroup(GG);
}


/******************************************************/
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( nrows(L) != nrows(S) )
  {
    ERROR("Incompatible matrices!");
  }

  def H = attrib(L, attrGroupHNF);
  if( typeof(H) != "intmat")
  {
    attrib(L, attrGroupHNF, hermiteNormalForm(L));
  } else { kill H; }

  def HH = attrib(L, attrGroupSNF);
  if( 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;
  // And now a quick-and-dirty fix of Singular inability to handle attribs of attribs:
  // For the use of a group as an attribute for multigraded rings
  G[3] = attrib(L, attrGroupHNF); 
  G[4] = attrib(L, attrGroupSNF);
  

  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;

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

}

/******************************************************/
static proc areIsomorphicGroups(def G, def H)
"USAGE: areIsomorphicGroups(G, H); G and H are groups
PURPOSE: Check whether G and H define isomorphic groups.
RETURN: int, 1 for TRUE, 0 otherwise
EXAMPLE: example areIsomorphicGroups; shows an example
"
{
  ERROR("areIsomorphicGroups: Not yet implemented!");
  return (1); // TRUE
}
example
{
  "EXAMPLE:"; echo=2;

  // TODO!

}

/******************************************************/
proc isGroup(def G)
"USAGE: isGroup(G); G a list
PURPOSE: checks whether G is a valid group
NOTE: G should be created by createGroup
(or createQuotientGroup, createTorsionFreeGroup)
RETURN: int, 1 if G is a valid group and 0 otherwise
EXAMPLE: example isGroup; shows an example
"
{
  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( defined(a) ){ if(typeof(a) == "int") { return(a); } }


  if( (size(G) != 2) && (size(G) != 4) ){ 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);
}
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

  isGroup(G);
  
  printGroup(G);

}



/******************************************************/
proc setBaseMultigrading(intmat M, list #)
"USAGE: setBaseMultigrading(M[, G]); M is an integer 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.
RETURN: nothing
EXAMPLE: example setBaseMultigrading; shows an example
"
{
  string attrMgrad   = "mgrad";
  string attrGradingGroup = "gradingGroup";

  int i = 1;
  if( size(#) >= i )
  {
    def a = #[i];
    if( typeof(a) == "intmat" )
    {
      def L = createGroup(M, a);
      i++;
    }

    if( isGroup(a) )
    {
      def L = a;

      if( !isSublattice(M, L[1]) )
      {
        ERROR("Multigrading is not contained in the grading group!");
      }
      i++;
    }
    if( i == 1 ){ ERROR("Wrong arguments: no group given?"); }
    kill a;
  }
  else
  {
    def L = createTorsionFreeGroup(M);
  }


  attrib(basering, attrMgrad, M);
  attrib(basering, attrGradingGroup, L);

  ideal Q = ideal(basering); // quotient ideal is assumed to be a GB!
  if( !isHomogeneous(Q, "checkGens") ) // 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;

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

  // GradingGroup:
  intmat L[3][2] =
    1, 1,
    1, 3,
    1, 5;

  // attaches M & L to R (==basering):
  setBaseMultigrading(M, L); // Grading: Z^3/L

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

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

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

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

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

  // Torsion:
  intmat L[2][1] =
    0,
    2;

  // attaches M & L to R (==basering):
  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)

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

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

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

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

  // Torsion:
  intmat L[1][1] =
    0;

  // attaches M & L to R (==basering):
  setBaseMultigrading(M); // Grading: Z^3

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

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


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


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

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

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

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

  // Grading group:
  intmat L[3][2] =
    1, 1,
    1, 3,
    1, 5;

  // attaches M & L to R (==basering):
  setBaseMultigrading(M, L); // Grading: Z^3/L

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

  kill L, M;

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

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

  // Grading group:
  intmat L[2][1] =
    0,
    2;

  // attaches M & L to R (==basering):
  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)

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

  kill L, M;

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

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

  // Grading group:
  intmat L[1][1] =
    0;

  // attaches M & L to R (==basering):
  setBaseMultigrading(M); // Grading: Z^3

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


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;

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

  if( i == 1 )
  {
    def R = basering;
  }

  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;

  // 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

  def G = getGradingGroup();

  printGroup( G );

  G[1] == M; G[2] == L;

  kill L, M, G;

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

  def G = getGradingGroup();

  printGroup( G );

  G[1] == M; G[2] == L;

  kill L, M, G;
  // ----------- 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

  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 )
  {
    def a = #[i];
    if( ( typeof(a) == "ring" ) or ( typeof(a) == "qring" ) )
    {
      i++;
    }
    kill a;
  }

  string attrGradingGroupHNF = "hermite";
  string attrGradingGroupSNF = "smith";

  def G = getGradingGroup(#);

//  printGroup(G);



  def T = G[2]; 

  if( size(#) >= i )
  {
    def a = #[i];

    if( typeof(a) != "string" )
    {
      ERROR("Sorry wrong arguments!");
    }
    
    if( a == "hermite" )
    {
      def M = attrib(T, attrGradingGroupHNF);
      if( typeof(M) != "intmat" )
      {
        if( size(G) > 2 )
        {
          M = G[3];
        } else
        {
          M = hermiteNormalForm(T);
        }
      }
      return (M);
    }

    if( a == "smith" )
    {
      def M = attrib(T, attrGradingGroupSNF);
      if( typeof(M) != "intmat" )
      {
        if( size(G) > 2 )
        {
          M = G[4];
        } else
        {
          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(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)
"USAGE: getGradedGenerator(M, i), 'M' module/ideal, 'i' int
RETURN: returns the i-th generator of M, endowed with the module grading from M
EXAMPLE: example getGradedGenerator; shows an example
"
{
  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");
}
example
{
  "EXAMPLE:"; echo=2;

  ring r = 0,(x,y,z,w),dp;
  intmat MM[2][4]=
                  1,1,1,1,
                  0,1,3,4;
  setBaseMultigrading(MM);
  
  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);


  intmat v[2][nrows(M)]=
                        1,
                        0;

  M = setModuleGrading(M, v);

  isHomogeneous(M);
  "Multidegrees: "; print(multiDeg(M));

  // Let's compute syzygies!
  def S = multiDegSyzygy(M); S;
  "Module Units Multigrading: "; print( getModuleGrading(S) );
  "Multidegrees: "; print(multiDeg(S));

  isHomogeneous(S);

  // same as S[1] together with the induced module weighting
  def v = getGradedGenerator(S, 1);
  print(v);
  print(setModuleGrading(v));
  
  isHomogeneous(v);

  isHomogeneous(S[1]);
}

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

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

  def V = attrib(m, attrModuleGrading);

  if( typeof(V) != "intmat" )
  {
    if( (typeof(m) == "ideal") or (typeof(m) == "poly") )
    {
      intmat M = getVariableWeights();
      intmat VV[nrows(M)][1];
      return (VV);
    }

    ERROR("Sorry: vector or module need module-grading-matrix! See 'getModuleGrading'.");
  }

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

  if( ncols(V) < nrows(m) )
  {
    ERROR("Sorry wrong width of V: " + string(ncols(V)));
  }

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

   ring R = 0, (x,y), dp;
   intmat M[2][2]=
     1, 1,
     0, 2;
   intmat T[2][5]=
     1,  2,  3,  4, 0,
     0, 10, 20, 30, 1;

   setBaseMultigrading(M, T);

   ideal I = x, y, xy^5;
   isHomogeneous(I);

   intmat V = multiDeg(I); print(V);

   module S = syz(I); print(S);

   S = setModuleGrading(S, V);

   getModuleGrading(S) == V;

   vector v = getGradedGenerator(S, 1);
   getModuleGrading(v) == V;
   isHomogeneous(v);
   print( multiDeg(v) );

   isHomogeneous(S);
   print( multiDeg(S) );
}

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

  intmat R = getVariableWeights();

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

  attrib(m, attrModuleGrading, G);
  return(m);
}
example
{
  "EXAMPLE:"; echo=2;

   ring R = 0, (x,y), dp;
   intmat M[2][2]=
     1, 1,
     0, 2;
   intmat T[2][5]=
     1,  2,  3,  4, 0,
     0, 10, 20, 30, 1;

   setBaseMultigrading(M, T);

   ideal I = x, y, xy^5;
   intmat V = multiDeg(I);

   // V == M; modulo T
   print(V);

   module S = syz(I);

   S = setModuleGrading(S, V);
   getModuleGrading(S) == V;

   print(S);

   vector v = getGradedGenerator(S, 1);
   getModuleGrading(v) == V;

   print( multiDeg(v) );

   isHomogeneous(S);

   print( multiDeg(S) );
}


proc multiDegTensor(module m, module n)
"
USAGE: multiDegTensor(m, n), m,n  modules or matrices.
PURPOSE: Computes the multigraded tensor product of to multigraded modules.
RETURN: A module.
EXAMPLE: example multiDegTensor; shows an example
"
{
  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 = multiDegTensor(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 = multiDegTensor(mm, nn);
print(mmtnn);
getModuleGrading(mmtnn);
matrix a = mmtnn;
matrix b = tensorMod(mm, nn);
print(a-b);

}

proc multiDegTor(int i, module m, module n)
{
  def res = multiDegResolution(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 = multiDegTensor(freedim3, m);
    //"mimag ok.";
    module mf = multiDegTensor(l[i], freedim);
    //"mf ok.";
    module mim1 = multiDegTensor(freedim2 ,m);
    module mim2 = multiDegTensor(l[i+1],freedim);
    //"mim1+2 ok.";
    module mker = multiDegModulo(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(multiDegModulo(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 = multiDegSyzygy(i);
module rt = Tor(2,m,m);
module multiDegT = multiDegTor(2,m,m);
print(matrix(rt)-matrix(multiDegT));
/*
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");
multiDegTor(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 = multiDegTensor(mm, nn);
*/
}


/******************************************************/
proc isGroupHomomorphism(def L1, def L2, intmat A)
"USAGE: isGoupHomomorphism(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 free.
RETURN: boolean, the result of the test
EXAMPLE: example isTorsionFree; shows an example
"
{
  intmat H = smithNormalForm(getLattice()); // TODO: ?cache it?  //******

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

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

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

   ring R = 0,(x,y),dp;
   intmat M[2][2]=
     1,0,
     0,1;
   intmat T[2][5]=
     1, 2, 3, 4, 0,
     0,10,20,30, 1;

   setBaseMultigrading(M,T);

   // Is the resulting group  free?
   isTorsionFree();

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

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

   kill R, A, B;
}


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]) div 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();
    intvec buf;

    for (i= 1; i <= sizeA; i++)
    {
      buf = intvec (m[i , 1..columns]);
      result= result+ list (buf);

    }
    return(result);
  }
  else
  {
    if(typeof(A) == "intmat")
    {
      A= system ("LLL", A);
      return(A);
    }
    else
    {
      ERROR("Unrecognized type.");
    }
  }
}
example
{
  "EXAMPLE:"; echo=2;

  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 #)
"USAGE: smithNormalForm(A[,opt]); intmat A
PURPOSE: Computes the Smith Normal Form of A
RETURN: if no optional argument is given: intmat, the Smith Normal Form of A,
otherwise: a list of 3 integer matrices P, D Q, such that D == P*A*Q.
EXAMPLE: example smithNormalForm; shows an example
"
{
  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] div 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 div v[1];
      transform[k+1,k+1] = b div 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;

  print( smithNormalForm(A) );

  list l = smithNormalForm(A, 5);

  l;

  l[1]*A*l[3];

  det(l[1]);
  det(l[3]);
}


/******************************************************/
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 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(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(A[row,column] == 0)
        {
          row++;
          continue;
        }
      for(j = column+1; j<=cc; j++)
        {
          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] div gcdvec[1];
              transform[column, column] = gcdvec[2];
              transform[column,j] = -v2[row] div gcdvec[1];
              transform[j,column] = gcdvec[3];
              q = q*transform;
              A = A*transform;
              // A;
            }
        }
      if(A[row,column]<0)
        {
          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]) div A[row, column];
          if(A[row,j]<0){
            transform[column,j]=transform[column,j]+1;}
          q = q*transform;
          A = A*transform;
        }
      }
      column++;
    }
  if(size(#) > 0){
    return(list(A, q));
  }
  return(A);
}
example
{
  "EXAMPLE:"; echo=2;

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

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

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

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

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

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

proc areZeroElements(intmat m, list #)
"USAGE: areZeroElements(D, [T]); intmat D, group T
PURPOSE: For a integer matrix D, considered column-wise as a set of
   integer vecors representing the multidegree of some polynomial
   or vector this method checks whether all these multidegrees
   are contained in the grading group
   group (either set globally or given as an optional argument),
i.e. if they all are zero in the multigrading.
EXAMPLE: example areZeroElements; shows an example
"
{
  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
{
  "EXAMPLE:"; echo=2;

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

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

  intmat L[2][1]=
    2,
    2;

  setBaseMultigrading(S,L);

  poly a = 1;
  poly b = xyz;

  ideal I = a, b;
  print(multiDeg(I));
  
  intmat m[5][2]=multiDeg(a),multiDeg(b); m=transpose(m);

  print(multiDeg(a));
  print(multiDeg(b));

  print(m);
 
  areZeroElements(m);

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

  areZeroElements(m,LL);
}


/******************************************************/
proc isZeroElement(intvec mdeg, list #)
"USAGE: isZeroElement(d, [T]); intvec d, group T
PURPOSE: For a integer vector 'd' 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
"
{
  int i = 1;
  if( size(#) >= i )
  {
    def a = #[1];
    if( typeof(a) == "intmat" )
    {
      intmat H = hermiteNormalForm(a);
      i++;
    }
    if( typeof(a) == "list" )
    {
      list L = a;
      intmat H = attrib(L, "hermite"); // todo
      i++;
    }
    kill a;
  }
  
  if( i == 1 )
  {
    intmat H = getLattice("hermite");
  }

  int x, k, row;

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

  int rr = nrows(mdeg);
  row = 1;
  intvec v;
  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]) div v[row]*v;
      }
    }
  }
  return( mdeg == 0 );

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

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

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

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

  setBaseMultigrading(g,t);

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

  intvec v1 = multiDeg(a) - multiDeg(b);
  v1;
  isZeroElement(v1);
  isZeroElement(v1, tt);

  intvec v2 = multiDeg(a) - multiDeg(c);
  v2;
  isZeroElement(v2);
  isZeroElement(v2, tt);

  intvec v3 = multiDeg(e) - multiDeg(f);
  v3;
  isZeroElement(v3);
  isZeroElement(v3, tt);

  intvec v4 = multiDeg(c) - multiDeg(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
"
{
  int i = 1;
  if( size(#) >= i )
  {
    def a = #[1];
    if( typeof(a) == "intmat" )
    {
      intmat grad = a;
      i++;
    }
    kill a;
  }

  if( i == 1 )
  {
    intmat grad = getVariableWeights();
  }

  intmat newgg[nrows(grad)][size(f)-1];
  int 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 oldlat=getLattice(pre);

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

  intmat newgrad[pp][np];

  //This will set the finest grading on the image ring. We will proceed by coarsening this grading until f becomes homogeneous. 
  for(i=1;i<=np;i++){ newgrad[p+i,i]=1;}

  //newgrad;



  list newlat;
  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];}
  }
  
  // This will make the images of homogeneous elements homogeneous, namely the variables of the preimage ring.
  for(i=1;i<=n;i++){
    im=f[i];
    //im;
    toadd = defineHomogeneous(im, newgrad);
    newlat=insert(newlat,toadd);
    columns=columns+ncols(toadd);

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

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

  //If the image ring is a quotient ring by some ideal, we have to coarsen the grading in order to make the ideal homogeneous.
  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);
      newlat=insert(newlat,toadd);
    }
  }

  //The grading group of the preimage ring might not have been torsion free. We have to add this torsion to the grading group of the image ring.
  intmat imofoldlat[pp][ncols(oldlat)];
  for(i=1; i<=nrows(oldlat);i++){
    for(j=1; j<=ncols(oldlat); j++){
      imofoldlat[i,j]=oldlat[i,j];
    }
  }

  columns=columns+ncols(oldlat);
  newlat=insert(newlat, imofoldlat);

  intmat gragr[pp][columns];
  columns=0;
  for(k=1;k<=size(newlat);k++){
    for(i=1;i<=pp;i++){
      for(j=1;j<=ncols(newlat[k]);j++){gragr[i,j+columns]=newlat[k][i,j];}
    }
    columns=columns+ncols(newlat[k]);
  }
  
  //The following is just for reducing the size of the matrices.
  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 equalMultiDeg(intvec exp1, intvec exp2, list #)
"USAGE: equalMultiDeg(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 equalMultiDeg; 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;


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

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

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

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

  equalMultiDeg(leadexp(e), leadexp(f));

}



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

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

  return (T == Z); // no grading group!
}


/******************************************************/
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
"
{
  if( (typeof(a) == "poly") or (typeof(a) == "vector") )
  {
    return ( size(multiDegPartition(a)) <= 1 )
  }

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

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

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

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

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

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

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

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

  setBaseMultigrading(M,T);

  attrib(r);

  poly f = x-yz;

  multiDegPartition(f);
  print(multiDeg(_));

  isHomogeneous(f);   // f: is not homogeneous

  poly g = 1-xy2z3;
  isHomogeneous(g); // g: is homogeneous
  multiDegPartition(g);

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

  setBaseMultigrading(M,T);

  f;
  isHomogeneous(f);
  multiDegPartition(f);

  // ---------------------
  g;
  isHomogeneous(g);
  multiDegPartition(g);

  kill r, T, M;

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

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

  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:
  multiDegPartition(x2 - y3 -xy +z);
  multiDegPartition(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;
  multiDegPartition(v1);
  print( multiDeg(_) );

  v2 = setModuleGrading(N[2], V); v2;
  multiDegPartition(v2);
  print( multiDeg(_) );

  N = setModuleGrading(N, V);
  isHomogeneous(N);
  print( multiDeg(N) );

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

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

  v1 = setModuleGrading(N[1], V); v1;
  multiDegPartition(v1);
  print( multiDeg(_) );

  v2 = setModuleGrading(N[2], V); v2;
  multiDegPartition(v2);
  print( multiDeg(_) );

  N = setModuleGrading(N, V);
  isHomogeneous(N);
  print( multiDeg(N) );

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

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

  N = gen(1) + x * gen(2), z*gen(3);
  N = setModuleGrading(N, V); print(N);
  isHomogeneous(N);
  print( multiDeg(N) );
  v1 = getGradedGenerator(N,1); print(v1);
  multiDegPartition(v1);
  print( multiDeg(_) );
  N = setModuleGrading(N, V); print(N);
  isHomogeneous(N);
  print( multiDeg(N) );
}

/******************************************************/
proc multiDeg(def A)
"USAGE: multiDeg(A); polynomial/vector/ideal/module A
PURPOSE: compute multidegree
EXAMPLE: example multiDeg; shows an example
"
{
  def a = attrib(A, "grad");
  if( typeof(a) == "intvec" || typeof(a) == "intmat" )
  {
    return (a);
  }

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

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

  int r = nrows(M);

  if( (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 )
  {
    intvec v; v[r] = 0;
    return (v);
  }

  intvec m; m[r] = 0;

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

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

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

    return(m);
  }


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

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

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

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

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

      A = A - lead(A);
    }

    return(m);
  }

  int i, j; intvec d;

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

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

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

    for( i = ncols(A); i > 0; i-- )
    {
      v = getGradedGenerator(A, i);

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

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

    }

    return(G);
  }

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

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

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

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

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

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


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

  setBaseMultigrading(A, Ta);

  multiDeg( x*x*y );

  multiDeg( y*y*y*x );

  multiDeg( x*y + x + 1 );

  multiDegPartition(x*y + x + 1);

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

//  "ideal:";

  ideal I = y*x*x, x*y*y*y;
  print( multiDeg(I) );

  print ( multiDeg(ideal(0)) );
  print ( multiDeg(ideal(0,0,0)) );

//  "vectors:";

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

  multiDeg( setModuleGrading(y*y*y*gen(2), B ));
  multiDeg( setModuleGrading(x*x*gen(1), B ));


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

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

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

//  "module:";

  module D = x*gen(1), y*gen(2);
  D;
  D = setModuleGrading(D, B);
  print( multiDeg( D ) );


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

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

};





/******************************************************/
proc multiDegPartition(def p)
"USAGE: multiDegPartition(def p), p polynomial/vector
RETURNS: an ideal/module consisting of multigraded-homogeneous parts of p
EXAMPLE: example multiDegPartition; shows an example
"
{ // TODO: What about an ideal or module???

  if( typeof(p) == "poly" )
  {
    ideal I;
    poly mp, t, tt;
    intmat V;
  }
  else
  {
    if(  typeof(p) == "vector" )
    {
      module I;
      vector mp, t, tt;
      intmat V = getModuleGrading(p);
    }
    else
    {
      ERROR("Wrong ARGUMENT type!");
    }
  }

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

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

      while( size(tt) > 0 )
      {
        // TODO: we do not cache matrices (M,T,H,V), which remain the same :(
        // TODO: we need some low-level procedure with all these arguments...!
        if( equalMultiDeg( leadexp(tt), m, V  ) )
        {
          mp = mp + lead(tt); // "mp", mp;
        }
        else
        {
          t = t + lead(tt);  //  "t", t;
        }

        tt = tt - lead(tt);
      }

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

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

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

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

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

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

  setBaseMultigrading(g,t);

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

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

  multiDegPartition(v, B);
}



/******************************************************/
static proc unitMatrix(int n)
{
  intmat A[n][n];

  for( int i = n; i > 0; i-- )
  {
    A[i,i] = 1;
  }

  return (A);
}



/******************************************************/
static proc finestMDeg(def r)
"
USAGE: finestMDeg(r); ring r
RETURN: ring, r endowed with the finest multigrading
TODO: not yet...
"
{
  def save = basering;
  setring (r);

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

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

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



  if( n > 0)
  {

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

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

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

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

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

  //  ERROR("nope");

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

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

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

  finestMDeg(q);

}




/******************************************************/
static proc newMap(map F, intmat Q, list #)
"
USAGE: newMap(F, Q[, P]); map F, intmat Q[, intmat P]
PURPOSE: endowe the map F with the integer matrices P [and Q]
"
{
  attrib(F, "Q", Q);

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

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


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

  ring r = integer, (x), dp;


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

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

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

  if( @bf == 1 )
  {
    setring @save@;
  }

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

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



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

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

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

//----------------------------------------------------------------------
// calling 4ti2 and converting output
// Singular's string is too clumsy for this, hence we first prepare
// using standard unix commands
//----------------------------------------------------------------------


   j=system("sh","hilbert -q -n sing4ti2"); ////////// be quiet + no loggin!!!

   j=system("sh", "awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.hil " +
                "| sed s/[\\\ \\\t\\\v\\\f]/,/g " +
                "| sed s/,+/,/g|sed s/,,/,/g " +
                "| sed s/,,/,/g " +
                "> sing4ti2.converted" );


//----------------------------------------------------------------------
// reading output of 4ti2
//----------------------------------------------------------------------
   link ausg=":r sing4ti2.converted";
//--- last entry ideal(0) is used to tie the list to the basering
//--- it will not be processed any further

   string s = read(ausg);

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

   int k = 3;

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



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

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

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

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

  poly m = 1; int e;

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

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

  return (m);

}

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

  intmat M = getVariableWeights(R);

  //  print(M);

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

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

  intmat T = getLattice(R);

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

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

    nc = ncols(A);

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


    //      print ( AA );

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

    //      print(K);

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

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

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

    //      "!";      print(B);

  }


  //  print(A);



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

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

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

  ideal Q = (ideal(basering));

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

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

  attrib(I, "ZeroPart", J);

  return (I);

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

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

  intmat g1[2][2]=1,0,0,1;
  intmat 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);
  intvec v3=4:1;
  v3;
  multiDegBasis(v3);

  setBaseMultigrading(g1,l);
  multiDegBasis(v1);
  setBaseMultigrading(g2);
  multiDegBasis(v2);

  intmat M[2][2] = 1, -1, -1, 1;
  intvec d = -2, 2;

  setBaseMultigrading(M);

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

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

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

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

  intvec d = 4, 1;

  setBaseMultigrading(M, L);

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


  kill R, M, d;

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


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

  setBaseMultigrading(M); // BUG????

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



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


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

  intvec d = 4, 1;

  setBaseMultigrading(M, T); // BUG????

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


proc multiDegSyzygy(def I)
"USAGE: multiDegSyzygy(I); I is a ideal or a module
PURPOSE: computes the multigraded syzygy module of I
RETURNS: module, the syzygy module of I
NOTE: generators of I must be multigraded homogeneous
EXAMPLE: example multiDegSyzygy; shows an example
"
{
  if( isHomogeneous(I, "checkGens") == 0)
  {
    ERROR ("Sorry: inhomogeneous input!");
  }
  module S = syz(I);
  S = setModuleGrading(S, multiDeg(I));
  return (S);
}
example
{
  "EXAMPLE:"; echo=2;

  ring r = 0,(x,y,z,w),dp;
  intmat MM[2][4]=
    1,1,1,1,
    0,1,3,4;
  setBaseMultigrading(MM);
  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);


  intmat v[2][nrows(M)]=
    1,
    0;

  M = setModuleGrading(M, v);

  isHomogeneous(M);
  "Multidegrees: "; print(multiDeg(M));
  // Let's compute syzygies!
  def S = multiDegSyzygy(M); S;
  "Module Units Multigrading: "; print( getModuleGrading(S) );
  "Multidegrees: "; print(multiDeg(S));

  isHomogeneous(S);
}



proc multiDegModulo(def I, def J)
"USAGE: multiDegModulo(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 multiDegModulo; 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, multiDeg(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(multiDeg(h1));

  // Let's compute modulo(h1, h2):
  def K = multiDegModulo(h1, h2); K;

  "Module Units Multigrading: "; print( getModuleGrading(K) );
  "Multidegrees: "; print(multiDeg(K));

  isHomogeneous(K);
}


proc multiDegGroebner(def I)
"USAGE: multiDegGroebner(I); I is a poly/vector/ideal/module
PURPOSE: computes the multigraded standard/groebner basis of I
NOTE: I must be multigraded homogeneous
RETURNS: ideal/module, the computed basis
EXAMPLE: example multiDegGroebner; shows an example
"
{
  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));
  }

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

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

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

  setBaseMultigrading(MM);


  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);


  intmat v[2][nrows(M)]=
    1,
    0;

  M = setModuleGrading(M, v);


  /////////////////////////////////////////////////////////////////////////////
  // GB:
  M = multiDegGroebner(M); M;
  "Module Units Multigrading: "; print( getModuleGrading(M) );
  "Multidegrees: "; print(multiDeg(M));

  isHomogeneous(M);

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

  isHomogeneous(S);

  /////////////////////////////////////////////////////////////////////////////
  // GB:
  S = multiDegGroebner(S); S;
  "Module Units Multigrading: "; print( getModuleGrading(S) );
  "Multidegrees: "; print(multiDeg(S));

  isHomogeneous(S);
}


/******************************************************/
proc multiDegResolution(def I, int ll, list #)
"USAGE: multiDegResolution(I,l,[f]); I is poly/vector/ideal/module; l,f are integers
PURPOSE: computes the multigraded resolution of I of the length l,
or the whole resolution if l is zero. Returns minimal resolution if an optional
argument 1 is supplied
NOTE: input must have multigraded-homogeneous generators.
The returned list is truncated beginning with the first zero differential.
RETURNS: list, the computed resolution
EXAMPLE: example multiDegResolution; shows an example
"
{
  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], V);
  }

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

  return (L);


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

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

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

  setBaseMultigrading(M);


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

  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++)
  {
    multiDegPartition(A[j]);
  }

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

  m = setModuleGrading(m, v);

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

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

  S = multiDegGroebner(S); S;
  "Module Units Multigrading: "; print( getModuleGrading(S) );
  "Multidegrees: "; print(multiDeg(S));

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

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

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

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

  L = multiDegResolution(maxideal(1), 0, 1);

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

  kill v;


  def h = hilbertSeries(m);
  setring h;

  numerator1;
  factorize(numerator1);

  denominator1;
  factorize(denominator1);

  numerator2;
  factorize(numerator2);

  denominator2;
  factorize(denominator2);
}

/******************************************************/
proc hilbertSeries(def I)
"USAGE: hilbertSeries(I); I is poly/vector/ideal/module
PURPOSE: computes the multigraded Hilbert Series of I
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 mutually prime numerator2
and denominator2, quotients of which give the series.
EXAMPLE: example hilbertSeries; shows an example
"
{

  if( !isFreeRepresented() )
  {
    "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);

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

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

  int l = size(RES);

  list L; L[l + 1] = 0;

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

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

  l++;

  ring R = 0,(t_(1..n),s_(1..n)),dp;

  ideal units;
  for( i=n; i>=1; i--)
  {
    units[i] = (var(i) * var(n + i) - 1);
  }

  qring Q = std(units);

  // TODO: should not it be a quotient ring depending on Torsion???
  // I am not sure about what to do in the torsion case, but since
  // we want to evaluate the polynomial at certain points to get
  // a dimension we need uniqueness for this. I think we would lose
  // this uniqueness if switching to this torsion ring.

  poly monom, summand, numerator;
  poly denominator = 1;

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

      if(v >= 0)
      {
        monom = monom * (var(k)^(v));
      }
      else
      {
        monom = monom * (var(n+k)^(-v));
      }
    }

    if( monom == 1)
    {
      ERROR("Multigrading not positive.");
    }

    denominator = denominator * (1 - monom);
  }

  for( j = 1; j<= l; j++)
  {
    summand = 0;
    M = L[j];

    for( i = 1; i <= ncols(M); i++)
    {
      monom = 1;
      for( k = 1; k <= n; k++)
      {
        v = M[k,i];
        if( v > 0 )
        {
          monom = monom * (var(k)^v);
        }
        else
        {
          monom = monom * (var(n+k)^(-v));
        }
      }
      summand = summand + monom;
    }
    numerator = numerator - (-1)^j * summand;
  }

  if( denominator == 0 )
  {
    ERROR("Multigrading not positive.");
  }

  poly denominator1 = denominator;
  poly numerator1 = numerator;

  export denominator1;
  export numerator1;

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

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

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


  export denominator2;
  export numerator2;

  " ------------ ";
  "This proc returns a ring with polynomials called 'numerator1/2' and 'denominator1/2'!";
  "They represent the first and the second Hilbert Series.";
  "The s_(i)-variables are defined to be the inverse of the t_(i)-variables.";
  " ------------ ";

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

  ring r = 0,(x,y,z,w),dp;
  intmat g[2][4]=
    1,1,1,1,
    0,1,3,4;
  setBaseMultigrading(g);

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

  M = setModuleGrading(M, V);

  def h = hilbertSeries(M); setring h;

  factorize(numerator2);
  factorize(denominator2);

  kill g, h; setring r;

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

  setBaseMultigrading(g);

  ideal I = x^2, y, z^3;
  I = std(I);
  def L = multiDegResolution(I, 0, 1);

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

  multiDeg(I);
  def h = hilbertSeries(I); setring h;

  factorize(numerator2);
  factorize(denominator2);

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

  def h = hilbertSeries(I); setring h;

  factorize(numerator2);
  factorize(denominator2);

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

  def h = hilbertSeries(I); setring h;

  factorize(numerator2);
  factorize(denominator2);

  kill R, W, h;
  ////////////////////////////////////////////////
  // This is example 5.3.9. from Robbianos book.

  ring R = 0,(x,y,z,w),dp;
  intmat W[1][4] =
     1,1, 1,1;
  setBaseMultigrading(W);
  ideal I = z3,y3zw2,x2y4w2xyz2;

  hilb(std(I));

  def h = hilbertSeries(I); setring h;

  numerator1;
  denominator1;

  factorize(numerator2);
  factorize(denominator2);


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

  ideal I2 = x2,y2,z2; I2;

  hilb(std(I2));

  def h = hilbertSeries(I2); setring h;

  numerator1;
  denominator1;


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

  W = 2,2,2,2;

  setBaseMultigrading(W);

  getVariableWeights();

  intvec w = 2,2,2,2;

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

  kill w;


  def h = hilbertSeries(I2); setring h;


  numerator1; denominator1;
  kill h;


  kill R, W;

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

  ideal I;

  I = 1; I;

  hilb(std(I));

  def h = hilbertSeries(I); setring h;

  numerator1; denominator1;

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

  I = x; I;

  hilb(std(I));

  def h = hilbertSeries(I); setring h;

  numerator1; denominator1;

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

  I = x^5; I;

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

  def h = hilbertSeries(I); setring h;

  numerator1; denominator1;


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

  I = x^10; I;

  hilb(std(I));

  def h = hilbertSeries(I); setring h;

  numerator1; denominator1;

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

  module M = 1;

  M = setModuleGrading(M, W);


  hilb(std(M));

  def h = hilbertSeries(M); setring h;

  numerator1; denominator1;

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

  kill M; module M = x^5*gen(1);
//  intmat V[1][3] = 0; // TODO: this would lead to a wrong result!!!?
  intmat V[1][1] = 0; // all gen(i) of degree 0!

  M = setModuleGrading(M, V);

  hilb(std(M));

  def h = hilbertSeries(M); setring h;

  numerator1; denominator1;

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

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

  kill V;

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

  N = setModuleGrading(N, V);

  hilb(std(N));

  def h = hilbertSeries(N); setring h;

  numerator1; denominator1;

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


  module S = M + N;

  S = setModuleGrading(S, V);

  hilb(std(S));

  def h = hilbertSeries(S); setring h;

  numerator1; denominator1;

  kill h;

  kill V;
  kill R, W;

}

static proc evalHilbertSeries(def h, intvec v)
"
   TODO
   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 = multiDegBasis(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 B[1][2]=1,1;
  setBaseMultigrading(B);
  isPositive(B);
}

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

  example getVariableWeights;
  example getLattice;
  example getGradingGroup;
  example getModuleGrading;


  example multiDeg;
  example multiDegPartition;


  example hermiteNormalForm;
  example isHomogeneous;
  example isTorsionFree;
  example pushForward;
  example defineHomogeneous;

  example equalMultiDeg;
  example isZeroElement;

  example multiDegResolution;

  "// ******************* example hilbertSeries ************************//";
  example hilbertSeries;


// example multiDegBasis; // needs 4ti2!

  "The End!";
}


static proc multiDegTruncate(def M, intvec md)
{
  "d: ";
  print(md);

  "M: ";
  module LL = M; // + L for d+1
  LL;
  print(multiDeg(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 = multiDegBasis(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(multiDeg(L));

  L = multiDegModulo(L, LL);
  L = multiDegGroebner(L);
//  L = minbase(prune(L));

  "??????????";
  print(L);
  print(multiDeg(L));

  V = getModuleGrading(L);

  // take out other degrees
  for(i = ncols(L); i > 0; i-- )
  {
    if( !equalMultiDeg( multiDeg(getGradedGenerator(L, i)), md ) )
    {
      L[i] = 0;
    }
  }

  L = simplify(L, 2);
  L = setModuleGrading(L, V);
  print(L);
  print(multiDeg(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);
  multiDegTruncate(h1, multiDeg(x));
  multiDegTruncate(h1, multiDeg(y));

  "XY:";
  module h2 = ideal(x, y);
  h2 = setModuleGrading(h2, V);
  multiDegTruncate(h2, multiDeg(x));
  multiDegTruncate(h2, multiDeg(y));
  multiDegTruncate(h2, multiDeg(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 LL = latticeBasis(Mul);

  return(LL);
}
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
{
  "EXAMPLE:"; echo=2;

  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

  // another example
  intmat B[2][2] =
    1,2,
    1,2;

  int d = intRank(B);

  print(B);
  print(d); // Should be 1

  kill A, B, r, d;

}

/*****************************************************/

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: integer, 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 and g != -1 )
    {
      return(0);
    }

    kill B, v, R;

  }

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

  //ring R = 0,(x,y),dp;
  intmat S2[3][3]=
    0, 2, 3,
    0, 1, 1,
    3, 0, 2;

  intmat S1[3][2]=
    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(L); 
  print(B); // should result in a matrix whose columns generate the same lattice as  [1,2,3] and [0,3,6]:

   // another example
  intmat C[2][4] =
    1,1,2,0,
    2,3,4,0;

  // should result in a matrix whose
  // colums create the same  lattice as
  // [0,1],[1,0]
  intmat D = latticeBasis(C);

  print(D);

  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 whose columns
  // generate the same lattice as [-1,2,-1,0],[2,-3,0,1]
  intmat B = kernelLattice(LL);
  
  print(B);

  // another example
  intmat C[2][4] =
    1,0,2,0,
    0,1,2,0;

  // should result in a matrix whose
  // colums create the same  lattice as
  // [-2,-2,1,0], [0,0,0,1]
  intmat D = kernelLattice(C);

  print(D);

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

  intmat L = kernelLattice(Con);

  // 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 (3x2)-matrix with columns e.g. [1,1,-1] and [0,3,-2] (the generated lattice should be identical)
  print(preimageLattice(P,B));

  // another example
  intmat W[3][3] =
    1,0,0,
    0,1,1,
    0,2,0;

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

  // should be a (3x2)-matrix with columns e.g. [1,0,0] and [0,0,-1] (the generated lattice should be identical)
  print(preimageLattice(W,Z));

}

/******************************************************/
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);
  print(b);

  // another example

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

  // should be 1
  int c = isPrimitiveSublattice(B);
  print(c);

  kill A, b, B, c;
}

/******************************************************/
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(P) == 1 )
  {
    return(1);
  }

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

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

  // should be 0
  int b = isIntegralSurjective(A);
  print(b);

  // another example
  intmat B[2][3] =
    1,1,5,
    2,3,6;

  // should be 1
  int c = isIntegralSurjective(B);
  print(c);

  kill A, b, B, c;
}

/******************************************************/
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[1][n]; // U now is the n-dim zero-vector
    }
    else
    {
      // we want a matrix with column operations so we transpose
      intmat BB = transpose(B);
      list L = hermiteNormalForm(BB, "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 such that the corresponding lattice is created by
  // [-1, 2], [-2, 3], [-1, 0] and [0, 1]
  print(projectLattice(B));

  // another example

  intmat BB[4][2] =
    1,0,
    0,1,
    0,0,
    0,0;

  // should result in a (2x4)-matrix such that the corresponding lattice is created by
  // [0,0],[0,0],[1,0],[0,1]
  print(projectLattice(BB));

  // another example

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

  // should result in the (1x3)-matrix that consists of just zeros
  print(projectLattice(BBB));

}

/******************************************************/
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 (3x2)-matrix with columns
  //  e.g. [0, 3, 6], [-3, 0, 3] (the lattice should be the same)
  print(intersectLattices(A,B));

  // another example
  intmat C[2][3] =
    1,0,0,
    3,2,5;

  intmat D[2][3] =
    4,5,0,
    0,5,0;

  // should result in a (3x2)-matrix whose columns generate the
  // same lattice as [1,5], [0, 20]
  print(intersectLattices(C,D));
}

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

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 * det(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);

  // another example
  intmat C[2][2] =
    2,1,
    3,2;
  
  intmat D  = intInverse(C);

  // should be the unit matrix
  print(C * D);

  kill A, B, C, D;
}


/******************************************************/
static 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 group 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: primitiveSpan(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 primitiveSpan; 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[3][2] =
    1,4,
    2,5,
    3,6;

  // should return a (3x2)-matrix whose columns
  // generate the same lattice as [1, 2, 3] and [0, 1, 2] 
  intmat R = primitiveSpan(V);
  print(R);

  // another example
  intmat W[2][2] =
    1,0,
    0,1;

  // should return a (2x2)-matrix whose columns
  // generate the same lattice as [1, 0] and [0, 1] 
  intmat S = primitiveSpan(W);
  print(S);

  kill V, R, S, W;
}

/***********************************************************/