source: git/Singular/LIB/fpadim.lib @ 1bf4c54

///////////////////////////////////////////////////////
version="$Id$";
category="Noncommutative";
info="
LIBRARY: fpadim.lib     Algorithms for quotient algebras in the letterplace case
AUTHORS: Grischa Studzinski,       grischa.studzinski@rwth-aachen.de
@*       Viktor Levandovskyy,      viktor.levandovskyy@math.rwth-aachen.de
@*       Karim Abou Zeid,          karim.abou.zeid@rwth-aachen.de

Support: Joint projects LE 2697/2-1 and KR 1907/3-1 of the Priority Programme SPP 1489:
@* 'Algorithmische und Experimentelle Methoden in Algebra, Geometrie und Zahlentheorie'
@* of the German DFG

OVERVIEW: Given the free algebra A = K<x_1,...,x_n> and a (finite) Groebner basis
@*      GB = {g_1,..,g_w}, one is interested in the K-dimension and in the
@*      explicit K-basis of A/<GB>.
@*      Therefore one is interested in the following data:
@*      - the Ufnarovskij graph induced by GB
@*      - the mistletoes of A/<GB>
@*      - the K-dimension of A/<GB>
@*      - the Hilbert series of A/<GB>
@*
@*      The Ufnarovskij graph is used to determine whether A/<GB> has finite
@*      K-dimension. One has to check if the graph contains cycles.
@*      For the whole theory we refer to [ufna]. Given a
@*      reduced set of monomials GB one can define the basis tree, whose vertex
@*      set V consists of all normal monomials w.r.t. GB. For every two
@*      monomials m_1, m_2 in V there is a direct edge from m_1 to m_2, if and
@*      only if there exists x_k in {x_1,..,x_n}, such that m_1*x_k = m_2. The
@*      set M = {m in V | there is no edge from m to another monomial in V} is
@*      called the set of mistletoes. As one can easily see it consists of
@*      the endpoints of the graph. Since there is a unique path to every
@*      monomial in V the whole graph can be described only from the knowledge
@*      of the mistletoes. Note that V corresponds to a basis of A/<GB>, so
@*      knowing the mistletoes we know a K-basis. The name mistletoes was given
@*      to those points because of these miraculous value and the algorithm is
@*      named sickle, because a sickle is the tool to harvest mistletoes.
@*      For more details see [studzins]. This package uses the Letterplace
@*      format introduced by [lls]. The algebra can either be represented as a
@*      Letterplace ring or via integer vectors: Every variable will only be
@*      represented by its number, so variable one is represented as 1,
@*      variable two as 2 and so on. The monomial x_1*x_3*x_2 for example will
@*      be stored as (1,3,2). Multiplication is concatenation. Note that there
@*      is no algorithm for computing the normal form yet, but for our case it
@*      is not needed. Note that the name fpadim.lib is short for dimensions of
@*      finite presented algebras.
@*

REFERENCES:

@*   [ufna] Ufnarovskij: Combinatorical and asymptotic methods in algebra, 1990
@*   [lls] Levandovskyy, La Scala: Letterplace ideals and non-commutative
Groebner bases, 2009
@*   [studzins] Studzinski: Dimension computations in non-commutative,
associative algebras, Diploma thesis, RWTH Aachen, 2010

ASSUMPTIONS:
@* - basering is always a Letterplace ring
@* - all intvecs correspond to Letterplace monomials
@* - if you specify a different degree bound d,
d <= attrib(basering,uptodeg) holds
@* In the procedures below, 'iv' stands for intvec representation
and 'lp' for the letterplace representation of monomials

PROCEDURES:

lpGkDim(G);                        computes the Gelfand Kirillov dimension of A/<G>
ivDHilbert(L,n[,d]);       computes the K-dimension and the Hilbert series
ivDHilbertSickle(L,n[,d]); computes mistletoes, K-dimension and Hilbert series
ivDimCheck(L,n);           checks if the K-dimension of A/<L> is infinite
lpGlDimBound(G);           computes upper bound of global dimension of A/<G>
ivHilbert(L,n[,d]);        computes the Hilbert series of A/<L> in intvec format
ivKDim(L,n[,d]);           computes the K-dimension of A/<L> in intvec format
ivMis2Base(M);             computes a K-basis of the factor algebra
ivMis2Dim(M);              computes the K-dimension of the factor algebra
ivOrdMisLex(M);            orders a list of intvecs lexicographically
ivSickle(L,n[,d]);         computes the mistletoes of A/<L> in intvec format
ivSickleHil(L,n[,d]);      computes the mistletoes and Hilbert series of A/<L>
ivSickleDim(L,n[,d]);      computes the mistletoes and the K-dimension of A/<L>
lpDHilbert(G[,d,n]);       computes the K-dimension and Hilbert series of A/<G>
lpDHilbertSickle(G[,d,n]); computes mistletoes, K-dimension and Hilbert series
lpHilbert(G[,d,n]);        computes the Hilbert series of A/<G> in lp format
lpDimCheck(G);             checks if the K-dimension of A/<G> is infinite
lpKDim(G[,d,n]);           computes the K-dimension of A/<G> in lp format
lpMis2Base(M);             computes a K-basis of the factor algebra
lpMis2Dim(M);              computes the K-dimension of the factor algebra
lpOrdMisLex(M);            orders an ideal of lp-monomials lexicographically
lpSickle(G[,d,n]);         computes the mistletoes of A/<G> in lp format
lpSickleHil(G[,d,n]);      computes the mistletoes and Hilbert series of A/<G>
lpSickleDim(G[,d,n]);      computes the mistletoes and the K-dimension of A/<G>
sickle(G[,m,d,h]);         can be used to access all lp main procedures


ivL2lpI(L);           transforms a list of intvecs into an ideal of lp monomials
iv2lp(I);             transforms an intvec into the corresponding monomial
iv2lpList(L);         transforms a list of intmats into an ideal of lp monomials
iv2lpMat(M);          transforms an intmat into an ideal of lp monomials
lp2iv(p);             transforms a polynomial into the corresponding intvec
lp2ivId(G);           transforms an ideal into the corresponding list of intmats
lpId2ivLi(G);         transforms a lp-ideal into the corresponding list of intvecs

SEE ALSO: freegb_lib
";

LIB "freegb.lib"; //for letterplace rings
LIB "general.lib";//for sorting mistletoes

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


//--------------- auxiliary procedures ------------------

static proc allVars(list L, intvec P, int n)
"USAGE: allVars(L,P,n); L a list of intmats, P an intvec, n an integer
RETURN: int, 0 if all variables are contained in the quotient algebra, 1 otherwise
"
{int i,j,r;
  intvec V;
  for (i = 1; i <= size(P); i++) {if (P[i] == 1){ j = i; break;}}
  V = L[j][1..nrows(L[j]),1];
  for (i = 1; i <= n; i++) {if (isInVec(i,V) == 0) {r = 1; break;}}
  if (r == 0) {return(1);}
  else {return(0);}
}

static proc checkAssumptions(int d, list L)
"PURPOSE: Checks, if all the Assumptions are holding
"
{if (typeof(attrib(basering,"isLetterplaceRing"))=="string") {ERROR("Basering is not a Letterplace ring!");}
  if (d > attrib(basering,"uptodeg")) {ERROR("Specified degree bound exceeds ring parameter!");}
  int i;
  for (i = 1; i <= size(L); i++)
  {if (entryViolation(L[i], attrib(basering,"lV")))
    {ERROR("Not allowed monomial/intvec found!");}
  }
  return();
}

static proc createStartMat(int d, int n)
"USAGE: createStartMat(d,n); d, n integers
RETURN: intmat
PURPOSE:Creating the intmat with all normal monomials in n variables and of degree d to start with
NOTE:   d has to be > 0
"
{intmat M[(n^d)][d];
  int i1,i2,i3,i4;
  for (i1 = 1; i1 <= d; i1++)  //Spalten
  {i2 = 1; //durchlaeuft Zeilen
    while (i2 <= (n^d))
    {for (i3 = 1; i3 <= n; i3++)
      {for (i4 = 1; i4 <= (n^(i1-1)); i4++)
        {M[i2,i1] = i3;
          i2 = i2 + 1;
        }
      }
    }
  }
  return(M);
}

static proc createStartMat1(int n, intmat M)
"USAGE: createStartMat1(n,M); n an integer, M an intmat
RETURN: intmat, with all variables except those in M
"
{int i;
  intvec V,Vt;
  V = M[(1..nrows(M)),1];
  for (i = 1; i <= size(V); i++) {if (isInVec(i,V) == 0) {Vt = Vt,i;}}
  if (Vt == 0) {intmat S; return(S);}
  else {Vt = Vt[2..size(Vt)]; intmat S [size(Vt)][1]; S[1..size(Vt),1] = Vt; return(S);}
}

static proc entryViolation(intmat M, int n)
"PURPOSE:checks, if all entries in M are variable-related
"
{int i,j;
  for (i = 1; i <= nrows(M); i++)
  {for (j = 1; j <= ncols(M); j++)
    {if(!((1<=M[i,j])&&(M[i,j]<=n))) {return(1);}}
  }
  return(0);
}

static proc findDimen(intvec V, int n, list L, intvec P, list #)
"USAGE: findDimen(V,n,L,P,degbound); V,P intvecs, n, an integer, L a list,
@*      degbound an optional integer
RETURN: int
PURPOSE:Computing the K-dimension of the quotient algebra
"
{int degbound = 0;
  if (size(#) > 0) {if (#[1] > 0) {degbound = #[1];}}
  int dimen,i,j,w,it;
  intvec Vt,Vt2;
  module M;
  if (degbound == 0)
  {for (i = 1; i <= n; i++)
    {Vt = V, i; w = 0;
      for (j = 1; j<= size(P); j++)
      {if (P[j] <= size(Vt))
        {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
          if (isInMat(Vt2,L[j]) > 0) {w = 1; break;}
        }
      }
      if (w == 0)
      {vector Vtt;
        for (it = 1; it <= size(Vt); it++){Vtt = Vtt + Vt[it]*gen(it);}
        M = M,Vtt;
        kill Vtt;
      }
    }
    if (size(M) == 0) {return(0);}
    else
    {M = simplify(M,2);
      for (i = 1; i <= size(M); i++)
      {kill Vt; intvec Vt;
        for (j =1; j <= size(M[i]); j++){Vt[j] =  int(leadcoef(M[i][j]));}
        dimen = dimen + 1 + findDimen(Vt,n,L,P);
      }
      return(dimen);
    }
  }
  else
  {if (size(V) > degbound) {ERROR("monomial exceeds degreebound");}
    if (size(V) == degbound) {return(0);}
    for (i = 1; i <= n; i++)
    {Vt = V, i; w = 0;
      for (j = 1; j<= size(P); j++)
      {if (P[j] <= size(Vt))
        {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
          if (isInMat(Vt2,L[j]) > 0) {w = 1; break;}
        }
      }
      if (w == 0) {vector Vtt;
        for (it = 1; it <= size(Vt); it++){Vtt = Vtt + Vt[it]*gen(it);}
        M = M,Vtt;
        kill Vtt;
      }
    }
    if (size(M) == 0) {return(0);}
    else
    {M = simplify(M,2);
      for (i = 1; i <= size(M); i++)
      {kill Vt; intvec Vt;
        for (j =1; j <= size(M[i]); j++){Vt[j] =  int(leadcoef(M[i][j]));}
        dimen = dimen + 1 + findDimen(Vt,n,L,P,degbound);
      }
      return(dimen);
    }
  }
}

static proc findCycle(intvec V, list L, intvec P, int n, int ld, module M)
"USAGE:
RETURN: int, 1 if Ufn-graph contains a cycle, or 0 otherwise
PURPOSE:Searching the Ufnarovskij graph for cycles
"
{int i,j,w,r;intvec Vt,Vt2;
  int it, it2;
  if (size(V) < ld)
  {for (i = 1; i <= n; i++)
    {Vt = V,i; w = 0;
      for (j = 1; j <= size(P); j++)
      {if (P[j] <= size(Vt))
        {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
          if (isInMat(Vt2,L[j]) > 0)
          {w = 1; break;}
        }
      }
      if (w == 0) {r = findCycle(Vt,L,P,n,ld,M);}
      if (r == 1) {break;}
    }
    return(r);
  }
  else
  {j = size(M);
    if (j > 0)
    {
      intmat Mt[j][nrows(M)];
      for (it = 1; it <= j; it++)
      { for(it2 = 1; it2 <= nrows(M);it2++)
        {Mt[it,it2] = int(leadcoef(M[it2,it]));}
      }
      Vt = V[(size(V)-ld+1)..size(V)];
      //Mt; type(Mt);Vt;type(Vt);
      if (isInMat(Vt,Mt) > 0) {return(1);}
      else
      {vector Vtt;
        for (it =1; it <= size(Vt); it++)
        {Vtt = Vtt + Vt[it]*gen(it);}
        M = M,Vtt;
        kill Vtt;
        for (i = 1; i <= n; i++)
        {Vt = V,i; w = 0;
          for (j = 1; j <= size(P); j++)
          {if (P[j] <= size(Vt))
            {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
              //L[j]; type(L[j]);Vt2;type(Vt2);
              if (isInMat(Vt2,L[j]) > 0)
              {w = 1; break;}
            }
          }
          if (w == 0) {r = findCycle(Vt,L,P,n,ld,M);}
          if (r == 1) {break;}
        }
        return(r);
      }
    }
    else
    { Vt = V[(size(V)-ld+1)..size(V)];
      vector Vtt;
      for (it = 1; it <= size(Vt); it++)
      {Vtt = Vtt + Vt[it]*gen(it);}
      M = Vtt;
      kill Vtt;
      for (i = 1; i <= n; i++)
      {Vt = V,i; w = 0;
        for (j = 1; j <= size(P); j++)
        {if (P[j] <= size(Vt))
          {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
            //L[j]; type(L[j]);Vt2;type(Vt2);
            if (isInMat(Vt2,L[j]) > 0)
            {w = 1; break;}
          }
        }
        if (w == 0) {r = findCycle(Vt,L,P,n,ld,M);}
        if (r == 1) {break;}
      }
      return(r);
    }
  }
}


static proc findCycleDFS(int i, intmat T, intvec V)
"
PURPOSE:
this is a classical deep-first search for cycles contained in a graph given by an intmat
"
{
  intvec rV;
  int k,k1,t;
  int j = V[size(V)];
  if (T[j,i] > 0) {return(V);}
  else
  {
    for (k = 1; k <= ncols(T); k++)
    {
      t = 0;
      if (T[j,k] > 0)
      {
        for (k1 = 1; k1 <= size(V); k1++) {if (V[k1] == k) {t = 1; break;}}
        if (t == 0)
        {
          rV = V;
          rV[size(rV)+1] = k;
          rV = findCycleDFS(i,T,rV);
          if (rV[1] > -1) {return(rV);}
        }
      }
    }
  }
  return(intvec(-1));
}



static proc findHCoeff(intvec V,int n,list L,intvec P,intvec H,list #)
"USAGE: findHCoeff(V,n,L,P,H,degbound); L a list of intmats, degbound an integer
RETURN: intvec
PURPOSE:Computing the coefficient of the Hilbert series (upto degree degbound)
NOTE:   Starting with a part of the Hilbert series we change the coefficient
@*      depending on how many basis elements we found on the actual branch
"
{int degbound = 0;
  if (size(#) > 0){if (#[1] > 0){degbound = #[1];}}
  int i,w,j,it;
  int h1 = 0;
  intvec Vt,Vt2,H1;
  module M;
  if (degbound == 0)
  {for (i = 1; i <= n; i++)
    {Vt = V, i; w = 0;
      for (j = 1; j<= size(P); j++)
      {if (P[j] <= size(Vt))
        {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
          if (isInMat(Vt2,L[j]) > 0) {w = 1; break;}
        }
      }
      if (w == 0)
      {vector Vtt;
        for (it = 1; it <= size(Vt); it++){Vtt = Vtt + Vt[it]*gen(it);}
        M = M,Vtt;
        kill Vtt;
      }
    }
    if (size(M) == 0) {return(H);}
    else
    {M = simplify(M,2);
      for (i = 1; i <= size(M); i++)
      {kill Vt; intvec Vt;
        for (j =1; j <= size(M[i]); j++) {Vt[j] =  int(leadcoef(M[i][j]));}
        h1 = h1 + 1; H1 = findHCoeff(Vt,n,L,P,H1);
      }
      if (size(H1) < (size(V)+2)) {H1[(size(V)+2)] = h1;}
      else {H1[(size(V)+2)] = H1[(size(V)+2)] + h1;}
      H1 = H1 + H;
      return(H1);
    }
  }
  else
  {if (size(V) > degbound) {ERROR("monomial exceeds degreebound");}
    if (size(V) == degbound) {return(H);}
    for (i = 1; i <= n; i++)
    {Vt = V, i; w = 0;
      for (j = 1; j<= size(P); j++)
      {if (P[j] <= size(Vt))
        {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
          if (isInMat(Vt2,L[j]) > 0) {w = 1; break;}
        }
      }
      if (w == 0)
      {vector Vtt;
        for (it = 1; it <= size(Vt); it++){Vtt = Vtt + Vt[it]*gen(it);}
        M = M,Vtt;
        kill Vtt;
      }
    }
    if (size(M) == 0) {return(H);}
    else
    {M = simplify(M,2);
      for (i = 1; i <= size(M); i++)
      {kill Vt; intvec Vt;
        for (j =1; j <= size(M[i]); j++)
        {Vt[j] =  int(leadcoef(M[i][j]));}
        h1 = h1 + 1; H1 = findHCoeff(Vt,n,L,P,H1,degbound);
      }
      if (size(H1) < (size(V)+2)) { H1[(size(V)+2)] = h1;}
      else {H1[(size(V)+2)] = H1[(size(V)+2)] + h1;}
      H1 = H1 + H;
      return(H1);
    }
  }
}

static proc findHCoeffMis(intvec V, int n, list L, intvec P, list R,list #)
"USAGE: findHCoeffMis(V,n,L,P,R,degbound); degbound an optional integer, L a
@*      list of Intmats, R
RETURN: list
PURPOSE:Computing the coefficients of the Hilbert series and the Mistletoes all
@*      at once
"
{int degbound = 0;
  if (size(#) > 0) {if (#[1] > 0) {degbound = #[1];}}
  int i,w,j,h1;
  intvec Vt,Vt2,H1; int it;
  module M;
  if (degbound == 0)
  {for (i = 1; i <= n; i++)
    {Vt = V, i; w = 0;
      for (j = 1; j<= size(P); j++)
      {if (P[j] <= size(Vt))
        {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
          if (isInMat(Vt2,L[j]) > 0) {w = 1; break;}
        }
      }
      if (w == 0)
      {vector Vtt;
        for (it = 1; it <= size(Vt); it++){Vtt = Vtt + Vt[it]*gen(it);}
        M = M,Vtt;
        kill Vtt;
      }
    }
    if (size(M) == 0) {if (size(R) < 2){R[2] = list(V);} else {R[2] = R[2] + list(V);} return(R);}
    else
    {M = simplify(M,2);
      for (i = 1; i <= size(M); i++)
      {kill Vt; intvec Vt;
        for (j =1; j <= size(M[i]); j++)
        {Vt[j] =  int(leadcoef(M[i][j]));}
        if (size(R[1]) < (size(V)+2)) { R[1][(size(V)+2)] = 1;}
        else
        {R[1][(size(V)+2)] = R[1][(size(V)+2)] + 1;}
        R = findHCoeffMis(Vt,n,L,P,R);
      }
      return(R);
    }
  }
  else
  {if (size(V) > degbound) {ERROR("monomial exceeds degreebound");}
    if (size(V) == degbound)
    {if (size(R) < 2){R[2] = list (V);}
      else{R[2] = R[2] + list (V);}
      return(R);
    }
    for (i = 1; i <= n; i++)
    {Vt = V, i; w = 0;
      for (j = 1; j<= size(P); j++)
      {if (P[j] <= size(Vt))
        {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
          if (isInMat(Vt2,L[j]) > 0) {w = 1; break;}
        }
      }
      if (w == 0)
      {vector Vtt;
        for (it = 1; it <= size(Vt); it++){Vtt = Vtt + Vt[it]*gen(it);}
        M = M,Vtt;
        kill Vtt;
      }
    }
    if (size(M) == 0) {if (size(R) < 2){R[2] = list(V);} else {R[2] = R[2] + list(V);} return(R);}
    else
    {M = simplify(M,2);
      for (i = 1; i <= ncols(M); i++)
      {kill Vt; intvec Vt;
        for (j =1; j <= size(M[i]); j++)
        {Vt[j] =  int(leadcoef(M[i][j]));}
        if (size(R[1]) < (size(V)+2)) { R[1][(size(V)+2)] = 1;}
        else
        {R[1][(size(V)+2)] = R[1][(size(V)+2)] + 1;}
        R = findHCoeffMis(Vt,n,L,P,R,degbound);
      }
      return(R);
    }
  }
}


static proc findMisDim(intvec V,int n,list L,intvec P,list R,list #)
"USAGE:
RETURN: list
PURPOSE:Computing the K-dimension and the Mistletoes all at once
"
{int degbound = 0;
  if (size(#) > 0) {if (#[1] > 0) {degbound = #[1];}}
  int dimen,i,j,w;
  intvec Vt,Vt2; int it;
  module M;
  if (degbound == 0)
  {for (i = 1; i <= n; i++)
    {Vt = V, i; w = 0;
      for (j = 1; j<= size(P); j++)
      {if (P[j] <= size(Vt))
        {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
          if (isInMat(Vt2,L[j]) > 0) {w = 1; break;}
        }
      }
      if (w == 0)
      {vector Vtt;
        for (it = 1; it <= size(Vt); it++){Vtt = Vtt + Vt[it]*gen(it);}
        M = M,Vtt;
        kill Vtt;
      }
    }
    if (size(M) == 0)
    {if (size(R) < 2){R[2] = list (V);}
      else{R[2] = R[2] + list(V);}
      return(R);
    }
    else
    {M = simplify(M,2);
      for (i = 1; i <= size(M); i++)
      {kill Vt; intvec Vt;
        for (j =1; j <= size(M[i]); j++){Vt[j] =  int(leadcoef(M[i][j]));}
        R[1] = R[1] + 1; R = findMisDim(Vt,n,L,P,R);
      }
      return(R);
    }
  }
  else
  {if (size(V) > degbound) {ERROR("monomial exceeds degreebound");}
    if (size(V) == degbound)
    {if (size(R) < 2){R[2] = list (V);}
      else{R[2] = R[2] + list (V);}
      return(R);
    }
    for (i = 1; i <= n; i++)
    {Vt = V, i; w = 0;
      for (j = 1; j<= size(P); j++)
      {if (P[j] <= size(Vt))
        {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
          if (isInMat(Vt2,L[j]) > 0) {w = 1; break;}
        }
      }
      if (w == 0)
      {vector Vtt;
        for (it = 1; it <= size(Vt); it++){Vtt = Vtt + Vt[it]*gen(it);}
        M = M,Vtt;
        kill Vtt;
      }
    }
    if (size(M) == 0)
    {if (size(R) < 2){R[2] = list (V);}
      else{R[2] = R[2] + list(V);}
      return(R);
    }
    else
    {M = simplify(M,2);
      for (i = 1; i <= size(M); i++)
      {kill Vt; intvec Vt;
        for (j =1; j <= size(M[i]); j++){Vt[j] =  int(leadcoef(M[i][j]));}
        R[1] = R[1] + 1; R = findMisDim(Vt,n,L,P,R,degbound);
      }
      return(R);
    }
  }
}


static proc findmistletoes(intvec V, int n, list L, intvec P, list #)
"USAGE: findmistletoes(V,n,L,P,degbound); V a normal word, n the number of
@*      variables, L the GB, P the occuring degrees,
@*      and degbound the (optional) degreebound
RETURN:  list
PURPOSE:Computing mistletoes starting in V
NOTE:   V has to be normal w.r.t. L, it will not be checked for being so
"
{int degbound = 0;
  if (size(#) > 0) {if (#[1] > 0) {degbound = #[1];}}
  list R; intvec Vt,Vt2; int it;
  int i,j;
  module M;
  if (degbound == 0)
  {int w;
    for (i = 1; i <= n; i++)
    {Vt = V,i; w = 0;
      for (j = 1; j <= size(P); j++)
      {if (P[j] <= size(Vt))
        {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
          if (isInMat(Vt2,L[j]) > 0)
          {w = 1; break;}
        }
      }
      if (w == 0)
      {vector Vtt;
        for (it = 1; it <= size(Vt); it++){Vtt = Vtt + Vt[it]*gen(it);}
        M = M,Vtt;
        kill Vtt;
      }
    }
    if (size(M)==0) {R = V; return(R);}
    else
    {M = simplify(M,2);
      for (i = 1; i <= size(M); i++)
      {kill Vt; intvec Vt;
        for (j =1; j <= size(M[i]); j++){Vt[j] =  int(leadcoef(M[i][j]));}
        R = R + findmistletoes(Vt,n,L,P);
      }
      return(R);
    }
  }
  else
  {if (size(V) > degbound) {ERROR("monomial exceeds degreebound");}
    if (size(V) == degbound) {R = V; return(R);}
    int w;
    for (i = 1; i <= n; i++)
    {Vt = V,i; w = 0;
      for (j = 1; j <= size(P); j++)
      {if (P[j] <= size(Vt))
        {Vt2 = Vt[(size(Vt)-P[j]+1)..size(Vt)];
          if (isInMat(Vt2,L[j]) > 0){w = 1; break;}
        }
      }
      if (w == 0)
      {vector Vtt;
        for (it = 1; it <= size(Vt); it++){Vtt = Vtt + Vt[it]*gen(it);}
        M = M,Vtt;
        kill Vtt;
      }
    }
    if (size(M) == 0) {R = V; return(R);}
    else
    {M = simplify(M,2);
      for (i = 1; i <= ncols(M); i++)
      {kill Vt; intvec Vt;
        for (j =1; j <= size(M[i]); j++)
        {Vt[j] =  int(leadcoef(M[i][j]));}
        //Vt; typeof(Vt); size(Vt);
        R = R + findmistletoes(Vt,n,L,P,degbound);
      }
      return(R);
    }
  }
}

static proc growthAlg(intmat T, list #)
"
real algorithm for checking the growth of an algebra
"
{
  int s = 1;
  if (size(#) > 0) { s = #[1];}
  int j;
  int n = ncols(T);
  intvec NV,C; NV[n] = 0; int m,i;
  intmat T2[n][n] = T[1..n,1..n]; intmat N[n][n];
  if (T2 == N)
  {
    for (i = 1; i <= n; i++)
    {
      if (m <  T[n+1,i]) { m = T[n+1,i];}
    }
    return(m);
  }

  //first part: the diagonals
  for (i = s; i <= n; i++)
  {
    if (T[i,i] > 0)
    {
      if ((T[i,i] >= 1) && (T[n+1,i] > 0)) {return(-1);}
      if ((T[i,i] == 1) && (T[n+1,i] == 0))
      {
        T[i,i] = 0;
        T[n+1,i] = 1;
        return(growthAlg(T));
      }
    }
  }

  //second part: searching for the last but one vertices
  T2 = T2*T2;
  for (i = s; i <= n; i++)
  {
    if ((intvec(T[i,1..n]) <> intvec(0)) && (intvec(T2[i,1..n]) == intvec(0)))
    {
      for (j = 1; j <= n; j++)
      {
        if ((T[i,j] > 0) && (m < T[n+1,j])) {m = T[n+1,j];}
      }
      T[n+1,i] = T[n+1,i] + m;
      T[i,1..n] = NV;
      return(growthAlg(T));
    }
  }
  m = 0;

  //third part: searching for circles
  for (i = s; i <= n; i++)
  {
    T2 = T[1..n,1..n];
    C = findCycleDFS(i,T2, intvec(i));
    if (C[1] > 0)
    {
      for (j = 2; j <= size(C); j++)
      {
        T[i,1..n] = T[i,1..n] + T[C[j],1..n];
        T[C[j],1..n] = NV;
      }
      for (j = 2; j <= size(C); j++)
      {
        T[1..n,i] = T[1..n,i] + T[1..n,C[j]];
        T[1..n,C[j]] = NV;
      }
      T[i,i] = T[i,i] - size(C) + 1;
      m = 0;
      for (j = 1; j <= size(C); j++)
      {
        m = m + T[n+1,C[j]];
      }
      for (j = 1; j <= size(C); j++)
      {
        T[n+1,C[j]] = m;
      }
      return(growthAlg(T,i));
    }
    else {ERROR("No Cycle found, something seems wrong! Please contact the authors.");}
  }

  m = 0;
  for (i = 1; i <= n; i++)
  {
    if (m < T[n+1,i])
    {
      m = T[n+1,i];
    }
  }
  return(m);
}

static proc GlDimSuffix(intvec v, intvec g)
{
  //Computes the shortest r such that g is a suffix for vr
  //only valid for lex orderings?
  intvec r,gt,vt,lt,g2;
  int lg,lv,l,i,c,f;
  lg = size(g); lv = size(v);
  if (lg <= lv)
  {
    l = lv-lg;
  }
  else
  {
    l = 0; g2 = g[(lv+1)..lg];
    g = g[1..lv]; lg = size(g);
    c = 1;
  }
  while (l < lv)
  {
    vt = v[(l+1)..lv];
    gt = g[1..(lv-l)];
    lt = size(gt);
    for (i = 1; i <= lt; i++)
    {
      if (vt[i]<>gt[i]) {l++; break;}
    }
    if (lt <=i ) { f = 1; break;}
  }
  if (f == 0) {return(g);}
  r = g[(lv-l+1)..lg];
  if (c == 1) {r = r,g2;}
  return(r);
}

static proc isNormal(intvec V, list G)
{
  int i,j,k,l;
  k = 0;
  for (i = 1; i <= size(G); i++)
  {
    if ( size(G[i]) <= size(V) )
    {
      while ( size(G[i])+k <= size(V) )
      {
        if ( G[i] == V[(1+k)..size(V)] ) {return(1);}
      }
    }
  }
  return(0);
}

static proc findDChain(list L)
{
  list Li; int i,j;
  for (i = 1; i <= size(L); i++) {Li[i] = size(L[i]);}
  Li = sort(Li); Li = Li[1];
  return(Li[size(Li)]);
}

static proc isInList(intvec V, list L)
"USAGE: isInList(V,L); V an intvec, L a list of intvecs
RETURN: int
PURPOSE:Finding the position of V in L, returns 0, if V is not in M
"
{int i,n;
  n = 0;
  for (i = 1; i <= size(L); i++) {if (L[i] == V) {n = i; break;}}
  return(n);
}

static proc isInMat(intvec V, intmat M)
"USAGE: isInMat(V,M);V an intvec, M an intmat
RETURN: int
PURPOSE:Finding the position of V in M, returns 0, if V is not in M
"
{if (size(V) <> ncols(M)) {return(0);}
  int i;
  intvec Vt;
  for (i = 1; i <= nrows(M); i++)
  {Vt = M[i,1..ncols(M)];
    if ((V-Vt) == 0){return(i);}
  }
  return(0);
}

static proc isInVec(int v,intvec V)
"USAGE: isInVec(v,V); v an integer,V an intvec
RETURN: int
PURPOSE:Finding the position of v in V, returns 0, if v is not in V
"
{int i,n;
  n = 0;
  for (i = 1; i <= size(V); i++) {if (V[i] == v) {n = i; break;}}
  return(n);
}


static proc isPF(intvec P, intvec I)
"
PURPOSE:
checks, if a word P is a praefix of another word I
"
{
  int n = size(P);
  if (n <= 0 || P == 0) {return(1);}
  if (size(I) < n) {return(0);}
  intvec IP = I[1..n];
  if (IP == P) {return(1);}
  else {return(0);}
}

static proc isSF(intvec S, intvec I)
"
PURPOSE:
checks, if a word S is a suffix of another word I
"
{
  int n = size(S);
  if (n <= 0 || S == 0) {return(1);}
  int m = size(I);
  if (m < n) {return(0);}
  intvec IS = I[(m-n+1)..m];
  if (IS == S) {return(1);}
  else {return(0);}
}

static proc isIF(intvec IF, intvec I)
"
PURPOSE:
checks, if a word IF is an infix of another word I
"
{
  int n = size(IF);
  int m = size(I);

  if (n <= 0 || IF == 0) {return(1);}
  if (m < n) {return(0);}

  for (int i = 0; (n + i) <= m; i++){
    intvec IIF = I[(1 + i)..(n + i)];
    if (IIF == IF) {
      return(1);
    }
  }
  return(0);
}

proc ivL2lpI(list L)
"USAGE: ivL2lpI(L); L a list of intvecs
RETURN: ideal
PURPOSE:Transforming a list of intvecs into an ideal of Letterplace monomials
ASSUME: - Intvec corresponds to a Letterplace monomial
@*      - basering has to be a Letterplace ring
EXAMPLE: example ivL2lpI; shows examples
"
{checkAssumptions(0,L);
  int i; ideal G;
  poly p;
  for (i = 1; i <= size(L); i++)
  {p = iv2lp(L[i]);
    G[(size(G) + 1)] = p;
  }
  return(G);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  def R = makeLetterplaceRing(5);// constructs a Letterplace ring
  setring R; //sets basering to Letterplace ring
  intvec u = 1,1,2; intvec v = 2,1,3; intvec w = 3,1,1;
  // u = x^2y, v = yxz, w = zx^2 in intvec representation
  list L = u,v,w;
  ivL2lpI(L);// invokes the procedure, returns the ideal containing u,v,w
}

proc iv2lp(intvec I)
"USAGE: iv2lp(I); I an intvec
RETURN: poly
PURPOSE:Transforming an intvec into the corresponding Letterplace polynomial
ASSUME: - Intvec corresponds to a Letterplace monomial
@*      - basering has to be a Letterplace ring
NOTE:   - Assumptions will not be checked!
EXAMPLE: example iv2lp; shows examples
"
{if (I[1] == 0) {return(1);}
  int i = size(I);
  if (i > attrib(basering,"uptodeg")) {ERROR("polynomial exceeds degreebound");}
  int j; poly p = 1;
  for (j = 1; j <= i; j++) {if (I[j] > 0) { p = lpMult(p,var(I[j]));}} //ignore zeroes, because they correspond to 1
  return(p);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; //sets basering to Letterplace ring
  intvec u = 1,1,2; intvec v = 2,1,3; intvec w = 3,1,1;
  // u = x^2y, v = yxz, w = zx^2 in intvec representation
  iv2lp(u); // invokes the procedure and returns the corresponding poly
  iv2lp(v);
  iv2lp(w);
}

proc iv2lpList(list L)
"USAGE: iv2lpList(L); L a list of intmats
RETURN: ideal
PURPOSE:Converting a list of intmats into an ideal of corresponding monomials
ASSUME: - The rows of each intmat in L must correspond to a Letterplace monomial
@*      - basering has to be a Letterplace ring
EXAMPLE: example iv2lpList; shows examples
"
{checkAssumptions(0,L);
  ideal G;
  int i;
  for (i = 1; i <= size(L); i++){G = G + iv2lpMat(L[i]);}
  return(G);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  intmat u[3][1] = 1,1,2; intmat v[1][3] = 2,1,3; intmat w[2][3] = 3,1,1,2,3,1;
  // defines intmats of different size containing intvec representations of
  // monomials as rows
  list L = u,v,w;
  print(u); print(v); print(w); // shows the intmats contained in L
  iv2lpList(L); // returns the corresponding monomials as an ideal
}


proc iv2lpMat(intmat M)
"USAGE: iv2lpMat(M); M an intmat
RETURN: ideal
PURPOSE:Converting an intmat into an ideal of the corresponding monomials
ASSUME: - The rows of M must correspond to Letterplace monomials
@*      - basering has to be a Letterplace ring
EXAMPLE: example iv2lpMat; shows examples
"
{list L = M;
  checkAssumptions(0,L);
  kill L;
  ideal G; poly p;
  int i; intvec I;
  for (i = 1; i <= nrows(M); i++)
  { I = M[i,1..ncols(M)];
    p = iv2lp(I);
    G[size(G)+1] = p;
  }
  return(G);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  intmat u[3][1] = 1,1,2; intmat v[1][3] = 2,1,3; intmat w[2][3] = 3,1,1,2,3,1;
  // defines intmats of different size containing intvec representations of
  // monomials as rows
  iv2lpMat(u); // returns the monomials contained in u
  iv2lpMat(v); // returns the monomials contained in v
  iv2lpMat(w); // returns the monomials contained in w
}

proc lpId2ivLi(ideal G)
"USAGE: lpId2ivLi(G); G an ideal
RETURN: list
PURPOSE:Transforming an ideal into the corresponding list of intvecs
ASSUME: - basering has to be a Letterplace ring
EXAMPLE: example lpId2ivLi; shows examples
"
{
  int i,j,k;
  list M;
  checkAssumptions(0,M);
  for (i = 1; i <= size(G); i++) {M[i] = lp2iv(G[i]);}
  return(M);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal L = x(1)*x(2),y(1)*y(2),x(1)*y(2)*x(3);
  lpId2ivLi(L); // returns the corresponding intvecs as a list
}

proc lp2iv(poly p)
"USAGE: lp2iv(p); p a poly
RETURN: intvec
PURPOSE:Transforming a monomial into the corresponding intvec
ASSUME: - basering has to be a Letterplace ring
NOTE:   - Assumptions will not be checked!
EXAMPLE: example lp2iv; shows examples
"
{p = normalize(lead(p));
  intvec I;
  int i,j;
  if (deg(p) > attrib(basering,"uptodeg")) {ERROR("Monomial exceeds degreebound");}
  if (p == 1) {return(I);}
  if (p == 0) {ERROR("Monomial is not allowed to equal zero");}
  intvec lep = leadexp(p);
  for ( i = 1; i <= attrib(basering,"lV"); i++) {if (lep[i] == 1) {I = i; break;}}
  for (i = (attrib(basering,"lV")+1); i <= size(lep); i++)
  {if (lep[i] == 1)
    { j = (i mod attrib(basering,"lV"));
      if (j == 0) {I = I,attrib(basering,"lV");}
      else {I = I,j;}
    }
    else { if (lep[i] > 1) {ERROR("monomial has a not allowed multidegree");}}
  }
  if (I[1] == 0) {ERROR("monomial has a not allowed multidegree");}

  return(I);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  poly p = x(1)*x(2)*z(3); poly q = y(1)*y(2)*x(3)*x(4);
  poly w= z(1)*y(2)*x(3)*z(4)*z(5);
  // p,q,w are some polynomials we want to transform into their
  // intvec representation
  lp2iv(p); lp2iv(q); lp2iv(w);
}

proc lp2ivId(ideal G)
"USAGE: lp2ivId(G); G an ideal
RETURN: list
PURPOSE:Converting an ideal into an list of intmats,
@*      the corresponding intvecs forming the rows
ASSUME: - basering has to be a Letterplace ring
EXAMPLE: example lp2ivId; shows examples
"
{G = normalize(lead(G));
  intvec I; list L;
  checkAssumptions(0,L);
  int i,md;
  for (i = 1; i <= size(G); i++) { if (md <= deg(G[i])) {md = deg(G[i]);}}
  while (size(G) > 0)
  {ideal Gt;
    for (i = 1; i <= ncols(G); i++) {if (md == deg(G[i])) {Gt = Gt + G[i]; G[i] = 0;}}
    if (size(Gt) > 0)
    {G = simplify(G,2);
      intmat M [size(Gt)][md];
      for (i = 1; i <= size(Gt); i++) {M[i,1..md] = lp2iv(Gt[i]);}
      L = insert(L,M);
      kill M; kill Gt;
      md = md - 1;
    }
    else {kill Gt; md = md - 1;}
  }
  return(L);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  poly p = x(1)*x(2)*z(3); poly q = y(1)*y(2)*x(3)*x(4);
  poly w = z(1)*y(2)*x(3)*z(4);
  // p,q,w are some polynomials we want to transform into their
  // intvec representation
  ideal G = p,q,w;
  // define the ideal containing p,q and w
  lp2ivId(G); // and return the list of intmats for this ideal
}

// -----------------main procedures----------------------

proc lpNoetherian(ideal G) {
  // return 0 not noetherian
  // return 1 left
  // return 2 right
  // return 3 both

  G = lead(G);
  G = simplify(G, 2+4+8);

  // check special case 1
  int l = 0;
  for (int i = 1; i <= size(G); i++) {
    // find the max degree in G
    int d = deg(G[i]);
    if (d > l) {
      l = d;
    }

    // also if G is the whole ring return noetherian
    if (leadmonom(G[i]) == 1) {
      return(3);
    }
  }
  // if longest word has length 1 we handle it as a special case
  if (l == 1) {
    int n = attrib(basering, "lV"); // variable count
    int k = size(G);
    if (k == n) { // only the field left
      return(3); // every field is noetherian
    }
    if (k == n-1) { // V = {1} with loop
      return(3);
    }
    if (k <= n-2) { // V = {1} with more than one loop
      return(0);
    }
  }

  intmat UG = lpUfGraph(G);

  // check special case 2
  intmat zero[nrows(UG)][ncols(UG)];
  if (UG == zero) {
    return (3);
  }

  if (!imHasLoops(UG) && imIsUpRightTriangle(topologicalSort(UG))) {
    // UG is a DAG
    return (3);
  }

  // DFS from every vertex, if cycle is found, check every vertex for incomming/outcom
  intvec visited;
  visited[ncols(UG)] = 0;
  int inFlag, outFlag;
  for (int v = 1; v <= ncols(UG) && (inFlag + outFlag) != 3; v++) {
    int inOutFlags = inOrOutCommingEdgeInCycle(UG, v, visited, 0);
    if (inOutFlags == 1) {
      inFlag = 1;
    }
    if (inOutFlags == 2) {
      outFlag = 2;
    }
    if (inOutFlags == 3) {
      inFlag = 1;
      outFlag = 2;
    }
  }
  return (3 - inFlag - outFlag);
}

proc inOrOutCommingEdgeInCycle(intmat G, int v, intvec visited, intvec path) {
  // Mark the current vertex as visited
  visited[v] = 1;

  // Store the current vertex in path
  if (path[1] == 0) {
    path[1] = v;
  } else {
    path[size(path) + 1] = v;
  }

  int inFlag, outFlag;

  for (int w = 1; w <= ncols(G) && (inFlag + outFlag) != 3; w++) {
    if (G[v,w] == 1) {
      if (visited[w] == 1) {
        // new cycle
        if (v == w) {
          for (int u = 1; u <= ncols(G); u++) {
            if (G[v,u] && u != v) {
              outFlag = 2;
            }
            if (G[u,v] && u != v) {
              inFlag = 1;
            }
          }
        } else {
          for (int i = size(path); i >= 1; i--) { // for each vertex in the path
            // check for neighbors not directly next or prev in cycle
            for (int u = 1; u <= ncols(G); u++) {
              if (G[path[i],u] == 1) { // there is an edge to u
                if (path[i] != v) {
                  if (u != path[i+1]) { // and u is not the next element in the cycle
                    outFlag = 2;
                  }
                } else {
                  if (u != w) {
                    outFlag = 2;
                  }
                }
              }
              if (G[u,path[i]] == 1) { // there is an edge from u
                if (path[i] != w) {
                  if (u != path[i-1]) { // and u is not the previous element in the cylce
                    inFlag = 1;
                  }
                } else {
                  if (u != v) {
                    inFlag = 1;
                  }
                }
              }
            }
            if (path[i] == w) {
              break;
            }
          }
        }
      } else {
        int inOutFlags = inOrOutCommingEdgeInCycle(G, w, visited, path);
        if (inOutFlags == 1) {
          inFlag = 1;
        }
        if (inOutFlags == 2) {
          outFlag = 2;
        }
        if (inOutFlags == 3) {
          inFlag = 1;
          outFlag = 2;
        }
      }
    }
  }

  return (inFlag + outFlag);
}

proc lpIsSemiPrime(ideal G)
{
  G = lead(G);
  G = simplify(G, 2+4+8);

  // check special case 1
  int l = 0;
  for (int i = 1; i <= size(G); i++) {
    // find the max degree in G
    int d = deg(G[i]);
    if (d > l) {
      l = d;
    }

    // also if G is the whole ring
    if (leadmonom(G[i]) == 1) {
      return(1);
    }
  }
  // if longest word has length 1 we handle it as a special case
  if (l == 1) {
    return(1);
  }

  list VUG = lpUfGraph(G, 1);
  intmat UG = VUG[1]; // the Ufnarovskij graph
  ideal V = VUG[2]; // the vertices of UG (standard words with length = l-1)

  list LG = lpId2ivLi(G);
  list SW = ivStandardWordsUpToLength(LG, maxDeg(G));
  list LV = lpId2ivLi(V);

  // delete the 0 in SW
  int indexofzero = ivIndexOf(SW, 0);
  if (indexofzero > 0) { // should be always true when |SW| > 0
    SW = delete(SW, indexofzero);
  }

  // check if each monomial in SW is cyclic
  for (int i = 1; i <= size(SW); i++) {
    if (!isCyclicInUfGraph(UG, LV, SW[i])) {
      return (0);
    }
  }

  return (1);
}

// checks whether a monomial is a cyclic monomial
proc isCyclicInUfGraph(intmat UG, list LV, intvec u)
{
  if (ncols(UG) == 0) {return (0);} // UG is empty
  if (u == 0) {return (0);} // 0 is never cyclic

  int l = size(LV[1]) + 1;

  int s = size(u);
  if (s <= l - 1) {
    for (int i = 1; i <= size(LV); i++) {
      // for all vertices where u is a suffix
      if(isSF(u, LV[i])) {
        if (existsRoute(UG, i, i)) {
          return (1);
        }
      }
    }
  } else { // size(u) > l - 1
    int m = s - l + 1;

    // there must be a route from v0 to vm
    intvec v0 = u[1..(l-1)]; // first in route of u
    intvec vm = u[m+1..m+(l-1)]; // last in route of u

    int iv0 = ivIndexOf(LV, v0);
    int ivm = ivIndexOf(LV, vm);
    if (iv0 <= 0 || ivm <= 0) {
      ERROR("u is not a standard word");
    }

    return (existsRoute(UG, ivm, iv0));
  }

  return (0);
}

proc lpIsPrime(ideal G)
"USAGE: lpIsPrime(G); G an ideal in a Letterplace ring
RETURN: boolean
PURPOSE: Check whether R/<G> is prime, where R is the basering
ASSUME: - basering is a Letterplace ring
@*      - G is a Groebner basis
"
{
  G = lead(G);
  G = simplify(G, 2+4+8);

  // check special case 1
  int l = 0;
  for (int i = 1; i <= size(G); i++) {
    // find the max degree in G
    int d = deg(G[i]);
    if (d > l) {
      l = d;
    }

    // also if G is the whole ring
    if (leadmonom(G[i]) == 1) {
      return(1);
    }
  }
  // if longest word has length 1 we handle it as a special case
  if (l == 1) {
    return(1);
  }

  list VUG = lpUfGraph(G, 1);
  intmat UG = VUG[1]; // the Ufnarovskij graph
  ideal V = VUG[2]; // the vertices of UG (standard words with length = l-1)

  list LG = lpId2ivLi(G);
  list LV = lpId2ivLi(V);

  int n = ncols(UG);

  // 1) for each vi vj there exists a route from vi to vj (means UG is connected)
  for (int i = 1; i <= n; i++) {
    for (int j = 1; j <= n; j++) {
      if (!existsRoute(UG, i, j)) {
        return (0);
      }
    }
  }

  // 2) any standard word with length < l-1 is a suffix of a vertex
  list SW = ivStandardWordsUpToLength(LG, maxDeg(G) - 2); // < maxDeg - 1
  if (size(SW) > 0 && size(LV) == 0) {return (0);}
  for (int i = 1; i <= size(SW); i++) {
    // check if SW[i] is a suffix of some LV
    for (int j = 1; j <= size(LV); j++) {
      if (!isSF(SW[i], LV[j])) {
        if (j == size(LV)) {
          return (0);
        }
      } else {
        break;
      }
    }
  }

  return (1);
}
example {
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5);
  setring R;
  ideal G = x(1)*x(2), y(1)*y(2); // K<x,y>/<xx,yy> is prime
  lpIsPrime(G);
}

static proc existsRoute(intmat G, int v, int u, list #)
"USAGE: existsRoute(G,v,u); G a graph, v and u vertices
NOTE: don't pass anything to # (internal use for recursion)
@*    routes always have at least one edge
"
{
  int n = ncols(G);

  // init visited
  intvec visited;
  if (size(#) > 0) {
    if (v == u) {return (1);} // don't check on first call so |route| >= 1 holds
    visited = #[1];
  } else { // first call
    visited[n] = 0;
  }

  // mark current vertex as visited
  visited[v] = 1;

  // recursive DFS
  for (int i = 1; i <= n; i++) {
    if (G[v,i] && (!visited[i] || i == u)) { // i == u to allow routes from u to u
      if (existsRoute(G, i, u, visited)) {
        return (1);
      }
    }
  }

  return (0);
}

static proc lpUfGkDim(ideal G)
{
  G = lead(G);
  G = simplify(G, 2+4+8);

  // check special case 1
  int l = 0;
  for (int i = 1; i <= size(G); i++) {
    // find the max degree in G
    int d = deg(G[i]);
    if (d > l) {
      l = d;
    }

    // also if G is the whole ring return minus infinity
    if (leadmonom(G[i]) == 1) {
      return(-2); // minus infinity
    }
  }
  // if longest word has length 1 we handle it as a special case
  if (l == 1) {
    int n = attrib(basering, "lV"); // variable count
    int k = size(G);
    if (k == n) { // V = {1} no edges
      return(0);
    }
    if (k == n-1) { // V = {1} with loop
      return(1);
    }
    if (k <= n-2) { // V = {1} with more than one loop
      return(-1);
    }
  }

  int t = rtimer; // DEBUG
  intmat UG = lpUfGraph(G);
  printf("lpUfGraph took %p", rtimer - t); // DEBUG

  // check special case 2
  intmat zero[nrows(UG)][ncols(UG)];
  if (UG == zero) {
    return (0);
  }

  // check special case 3
  UG = topologicalSort(UG);

  if (imIsUpRightTriangle(UG)) {
    UG = eliminateZerosUpTriangle(UG);
    if (ncols(UG) == 0 || nrows(UG) == 0) { // when the diagonal was zero
      return (0)
    }
    return(UfGraphURTNZDGrowth(UG));
  }

  // otherwise count cycles in the Ufnarovskij Graph
  int t = rtimer; // DEBUG
  int gkdim = UfGraphGrowth(UG);
  printf("UfGraphGrowth took %p", rtimer - t); // DEBUG
  return(gkdim);
}

static proc UfGraphURTNZDGrowth(intmat UG) {
  // URTNZD = upper right triangle non zero diagonal
  for (int i = 1; i <= ncols(UG); i++) {
    UG[i,i] = 0; // remove all loops
  }
  intmat UGk = UG;
  intmat zero[nrows(UGk)][ncols(UGk)];
  int k = 1;
  while (UGk != zero) {
    UGk = UGk * UG;
    k++;
  }
  return (k);
}

proc imIsUpRightTriangle(intmat M) {
  for (int i = 1; i <= nrows(M); i++) {
    for (int j = 1; j < i; j++) {
      if(M[i,j] != 0) { return (0); }
    }
  }
  return (1);
}

static proc eliminateZerosUpTriangle(intmat G) {
  // G is expected to be an upper triangle matrix
  for (int i = ncols(G); i >= 1; i--) { // loop order is important because we delete entries
    if (G[i,i] == 0) { // i doesn't have a cycle
      for (int j = 1; j < i; j++) {
        if (G[j,i] == 1) { // j has an edge to i
          for (int k = i + 1; k <= nrows(G); k++) {
            if (G[i,k] == 1) {
              G[j,k] = G[i,k]; // give j all edges from i
            }
          }
        }
      }
      G = imDelRowCol(G,i,i); // remove vertex i
    }
  }
  return (G);
}

static proc imDelRowCol(intmat M, int row, int col) {
  // row and col are expected to be > 0
  int nr = nrows(M);
  int nc = ncols(M);
  intmat Mdel[nr - 1][nc - 1];
  for (int i = 1; i <= nr; i++) {
    for (int j = 1; j <= nc; j++) {
      if(i != row && j != col) {
        int newi = i;
        int newj = j;
        if (i > row) { newi = i - 1; }
        if (j > col) { newj = j - 1; }
        Mdel[newi,newj] = M[i,j];
      }
    }
  }
  return (Mdel);
}

static proc topologicalSort(intmat G) {
  // NOTE: ignores loops
  // NOTE: this takes O(|V^3|), can be optimized
  int n = ncols(G);
  for (int i = 1; i <= n; i++) { // only use the submat at i
    // find a vertex v in the submat at i with no incoming edges
    int v;
    for (int j = i; j <= n; j++) {
      int incoming = 0;
      for (int k = i; k <= n; k++) {
        if (k != j && G[k,j] == 1) {
          incoming = 1;
        }
      }
      if (incoming == 0) {
        v = j;
        break;
      } else {
        if (j == n) {
          // G contains at least one cycle, abort
          return (G);
        }
      }
    }

    // swap v and i
    if (v != i) {
      G = imPermcol(G, v, i);
      G = imPermrow(G, v, i);
    }
  }
  return (G);
}

static proc imPermcol (intmat A, int c1, int c2)
{
  intmat B = A;
  int k = nrows(B);
  B[1..k,c1] = A[1..k,c2];
  B[1..k,c2] = A[1..k,c1];
  return (B);
}

static proc imPermrow (intmat A, int r1, int r2)
{
  intmat B = A;
  int k = ncols(B);
  B[r1,1..k] = A[r2,1..k];
  B[r2,1..k] = A[r1,1..k];
  return (B);
}

static proc UfGraphGrowth(intmat UG)
{
  int n = ncols(UG); // number of vertices
  // iterate through all vertices

  intvec visited;
  visited[n] = 0;

  intvec cyclic;
  cyclic[n] = 0;

  int maxCycleCount = 0;
  for (int v = 1; v <= n; v++) {
    int cycleCount = countCycles(UG, v, visited, cyclic, 0);
    if (cycleCount == -1) {
      return(-1);
    }
    if (cycleCount > maxCycleCount) {
      maxCycleCount = cycleCount;
    }
  }
  return(maxCycleCount);
}

static proc countCycles(intmat G, int v, intvec visited, intvec cyclic, intvec path)
"USAGE: countCycles(G, v, visited, cyclic, path); G a Graph, v the vertex to
start. The parameter visited, cyclic and path should be 0.
RETURN: int
@*:     Maximal number of distinct cycles
PURPOSE: Calculate the maximal number of distinct cycles in a single path starting at v
ASSUME: Basering is a Letterplace ring
"
{
  // Mark the current vertex as visited
  visited[v] = 1;

  // Store the current vertex in path
  if (path[1] == 0) {
    path[1] = v;
  } else {
    path[size(path) + 1] = v;
  }

  int cycles = 0;
  for (int w = 1; w <= ncols(G); w++) {
    if (G[v,w] == 1) {
      if (visited[w] == 1) { // neuer zykel gefunden
        // 1. alle Knoten in path bis w ÃŒberprÃŒfen ob in cyclic
        for (int j = size(path); j >= 1; j--) {
          if(cyclic[path[j]] == 1) {
            // 1.1 falls ja return -1
            return (-1);
          }
          if (path[j] == w) {
            break;
          }
        }

        // 2. ansonsten cycles++
        for (int j = size(path); j >= 1; j--) {
          // 2.2 Kanten in diesem Zykel entfernen; Knoten cyclic
          if (j == size(path)) { // Sonderfall bei der ersten Iteration
            cyclic[v] = 1;
            G[v, w] = 0;
          } else {
            cyclic[path[j]] = 1;
            G[path[j], path[j+1]] = 0;
          }
          if (path[j] == w) {
            break;
          }
        }

        // 3. auf jedem dieser Knoten countCycles() aufrufen
        int maxCycleCount = 0;
        for (int j = size(path); j >= 1; j--) {
          int cycleCount = countCycles(G, path[j], visited, cyclic, path);
          if(cycleCount == -1) {
            return (-1);
          }
          if (cycleCount > maxCycleCount) {
            maxCycleCount = cycleCount;
          }
          if (path[j] == w) {
            break;
          }
        }
        if (maxCycleCount >= cycles) {
          cycles = maxCycleCount + 1;
        }
      } else {
        int cycleCount = countCycles(G, w, visited, cyclic, path);
        if (cycleCount == -1) {
          return(-1);
        }
        if (cycleCount > cycles) {
          cycles = cycleCount;
        }
      }
    }
  }
  // printf("Path: %s countCycles: %s", path, cycles); // DEBUG
  return(cycles);
}

static proc lpUfGraph(ideal G, list #)
"USAGE: lpUfGraph(G); G a set of monomials in a letterplace ring, # an optional parameter to return the vertex list when set
RETURN: intmat
PURPOSE: Constructs the Ufnarovskij graph induced by G
@*      the adjacency matrix of the Ufnarovskij graph induced by G
ASSUME: - basering is a Letterplace ring
@*      - G are the leading monomials of a Groebner basis
"
{
  int l = maxDeg(G);
  list LG = lpId2ivLi(G);
  list SW = ivStandardWords(LG, l - 1); // vertices
  int n = size(SW);
  intmat UG[n][n]; // Ufnarovskij graph
  for (int i = 1; i <= n; i++) {
    for (int j = 1; j <= n; j++) {
      // [Studzinski page 76]
      intvec v = SW[i];
      intvec w = SW[j];
      intvec v_overlap;
      intvec w_overlap;
      //TODO how should the graph look like when l - 1 = 0 ?
      if (l - 1 == 0) {
        ERROR("Ufnarovskij graph not implemented for l = 1");
      }
      if (l - 1 > 1) {
        v_overlap = v[2 .. l-1];
        w_overlap = w[1 .. l-2];
      }
      intvec vw = v;
      vw[l] = w[l-1];
      if (v_overlap == w_overlap && !ivdivides(LG, vw)) {
        UG[i,j] = 1;
      }
    }
  }
  if (size(#) > 0) {
    if (typeof(#[1]) == "int") {
      if (#[1] == 1) {
        list ret = UG;
        ret[2] = ivL2lpI(SW); // the vertices
        return (ret);
      }
    }
  }
  return (UG);
}

static proc maxDeg(ideal G)
{
  int l = 0;
  for (int i = 1; i <= size(G); i++) { // find the max degree in G
    int d = deg(G[i]);
    if (d > l) {
      l = d;
    }
  }
  return (l);
}

static proc ivStandardWords(list G, int length)
"ASSUME: G is simplified
"
{
  if (length <= 0) {
    list words;
    if (length == 0 && !ivdivides(G,0)) {
      words[1] = 0; // iv = 0 means monom = 1
    }
    return (words); // no standard words
  }
  int lV = attrib(basering, "lV"); // variable count
  list prevWords = ivStandardWords(G, length - 1);
  list words;
  for (int i = 1; i <= lV; i++) {
    for (int j = 1; j <= size(prevWords); j++) {
      intvec word = prevWords[j];
      word[length] = i;
      // assumes that G is simplified!
      if (!ivdivides(G, word)) {
        words = insert(words, word);
      }
    }
  }
  return (words);
}

static proc ivStandardWordsUpToLength(list G, int length)
"ASSUME: G is simplified
"
{
  list words = ivStandardWords(G,0);
  if (size(words) == 0) {return (words)}
  for (int i = 1; i <= length; i++) {
    words = words + ivStandardWords(G, i);
  }
  return (words);
}

static proc ivdivides(list G, intvec iv) {
  for (int k = 1; k <= size(G); k++) {
    if (isIF(G[k], iv)) {
      return (1);
    } else {
      if (k == size(G)) {
        return (0);
      }
    }
  }
}

static proc lpGraphOfNormalWords(ideal G)
"USAGE: lpGraphOfNormalWords(G); G a set of monomials in a letterplace ring
RETURN: intmat
PURPOSE: Constructs the graph of normal words induced by G
@*:      the adjacency matrix of the graph of normal words induced by G
ASSUME: - basering is a Letterplace ring
- G are the leading monomials of a Groebner basis
"
{
  // construct the Graph of normal words [Studzinski page 78]
  // construct set of vertices
  int v = attrib(basering,"lV"); int d = attrib(basering,"uptodeg");
  ideal V; poly p,q,w;
  ideal LG = lead(G);
  int i,j,k,b; intvec E,Et;
  for (i = 1; i <= v; i++){V = V, var(i);}
  for (i = 1; i <= size(LG); i++)
  {
    E = leadexp(LG[i]);
    if (E == intvec(0)) {V = V,monomial(intvec(0));}
    else
    {
      for (j = 1; j < d; j++)
      {
        Et = E[(j*v+1)..(d*v)];
        if (Et == intvec(0)) {break;}
        else {V = V, monomial(Et);}
      }
    }
  }
  V = simplify(V,2+4);
  printf("V = %p", V);


  // construct incidence matrix

  list LV = lpId2ivLi(V);
  intvec Ip,Iw;
  int n = size(V);
  intmat T[n+1][n];
  for (i = 1; i <= n; i++)
  {
    // printf("for1 (i=%p, n=%p)", i, n);
    p = V[i]; Ip = lp2iv(p);
    for (j = 1; j <= n; j++)
    {
      // printf("for2 (j=%p, n=%p)", j, n);
      k = 1; b = 1;
      q = V[j];
      w = lpNF(lpMult(p,q),LG);
      if (w <> 0)
      {
        Iw = lp2iv(w);
        while (k <= n)
        {
          // printf("while (k=%p, n=%p)", k, n);
          if (isPF(LV[k],Iw) > 0)
          {if (isPF(LV[k],Ip) == 0) {b = 0; k = n+1;} else {k++;}
          }
          else {k++;}
        }
        T[i,j] = b;
        //  print("Incidence Matrix:");
        // print(T);
      }
    }
  }
  return(T);
}

static proc lpNorGkDim(ideal G)
"USAGE: lpNorGkDim(G); G an ideal in a letterplace ring
RETURN: int
PURPOSE: Determines the Gelfand Kirillov dimension of A/<G>
@*:     -1 means it is infinite
ASSUME: - basering is a Letterplace ring
- G is a Groebner basis
"
{
  return(growthAlg(lpGraphOfNormalWords(G)));
}

proc lpGkDim(ideal G, list#)
"USAGE: lpGkDim(G); G an ideal in a letterplace ring, method an
@*       optional integer. method = 0 uses the Ufnarovskij Graph
@*       and method = 1 uses the Graph of normal words to determine
@*       the Gelfand Kirillov dimension
RETURN: int
PURPOSE: Determines the Gelfand Kirillov dimension of A/<G>
@*      -1 means it is positive infinite
@*      -2 means it is negative infinite
ASSUME: - basering is a Letterplace ring
- G is a Groebner basis
"
{
  int method = 0;
  if (size(#) > 0) {
    if (typeof(#[1])=="int") {
      method = #[1];
    }
  }

  if (method == 0) {
    return (lpUfGkDim(G));
  } else {
    return (lpNorGkDim(G));
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  R;
  setring R; // sets basering to Letterplace ring
  ideal I = z(1);//an example of infinite GK dimension
  lpGkDim(I);
  I = x(1),y(1),z(1); // gkDim = 0
  lpGkDim(I);
  I = x(1)*y(2), x(1)*z(2), z(1)*y(2), z(1)*z(2);//gkDim = 2
  lpGkDim(I);
}

proc lpGlDimBound(ideal G)
"USAGE: lpGlDimBound(I); I an ideal
RETURN: int, an upper bound for the global dimension, -1 means infinity
PURPOSE: computing an upper bound for the global dimension
ASSUME: - basering is a Letterplace ring, G is a reduced Gröbner Basis
EXAMPLE: example lpGlDimBound; shows example
NOTE: if I = LM(I), then the global dimension is equal the Gelfand
@*    Kirillov dimension if it is finite
@*    Global dimension should be 0 for A/G = K and 1 for A/G = K<x1...xn>
"
{
  G = simplify(G,2); // remove zero generators
  // NOTE: Gl should be 0 for A/G = K and 1 for A/G = K<x1...xn>
  // G1 contains generators with single variable in LM
  ideal G1;
  for (int i = 1; i <= size(G); i++) {
    if (ord(G[i]) < 2) { // single variable in LM
      G1 = insertGenerator(G1,G[i]);
    }
  }
  G1 = simplify(G1,2); // remove zero generators

  // G = NF(G,G1)
  for (int i = 1; i <= ncols(G); i++) { // do not use size() here
    G[i] = lpNF(G[i],G1);
  }
  G = simplify(G,2); // remove zero generators

  // delete variables in LM(G1) from the ring
  def save = basering;
  ring R = basering;
  if (size(G1) > 0) {
    while (size(G1) > 0) {
      if (attrib(R, "lV") > 1) {
        ring R = lpDelVar(lp2iv(G1[1])[1]);
        ideal G1 = imap(save,G1);
        G1 = simplify(G1, 2); // remove zero generators
      } else {
        // only the field is left (no variables)
        return(0);
      }
    }
    ideal G = imap(save, G); // put this here, because when save == R this call would make G = 0
  }

  // Li p. 184 if G = LM(G), then I = LM(I) and thus glDim = gkDim if it's finite
  for (int i = 1; i <= size(G); i++) {
    if (G[i] != lead(G[i])) {
      break;
    } else {
      if (i == size(G)) { // if last iteration
        print("Using gk dim"); // DEBUG
        int gkDim = lpGkDim(G);
        if (gkDim >= 0) {
          return (gkDim);
        }
      }
    }
  }

  intmat GNC = lpGraphOfNChains(G);

  // assuming GNC is connected

  // TODO: maybe loop+cycle checking could be done more efficiently?
  if (!imHasLoops(GNC) && imIsUpRightTriangle(topologicalSort(GNC))) {
    // GNC is a DAG
    intmat GNCk = GNC;
    intmat zero[1][ncols(GNCk)];
    int k = 1;
    // while first row isn't empty
    while (GNCk[1,1..(ncols(GNCk))] != zero[1,1..(ncols(zero))]) {
      GNCk = GNCk * GNC;
      k++;
    }
    // k-1 = number of edges in longest path starting from 1
    return (k-1);
  } else {
    // GNC contains loops/cycles => there is always an n-chain
    return (-1); // infinity
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal G = x(1)*x(2), y(1)*y(2),x(1)*y(2)*x(3); // ideal G contains a
  //Groebner basis
  lpGlDimBound(G); // invokes procedure with Groebner basis G
}

static proc imHasLoops(intmat A) {
  int n = ncols(A);
  for (int i = 1; i < n; i++) {
    if (A[i,i] == 1) {
      return (1);
    }
  }
  return (0);
}

static proc lpGraphOfNChains(ideal G) // G must be reduced
{
  list LG = lpId2ivLi(lead(G));
  int n = attrib(basering, "lV");
  int degbound = attrib(basering, "uptodeg");

  list V;
  for (int i = 0; i <= n; i++) {
    V[i+1] = i; // add 1 and all variables
  }
  for (int i = 1; i <= size(LG); i++) {
    intvec u = LG[i];
    for (int j = 2; j <= size(u); j++) {
      intvec v = u[j..size(u)];
      if (!contains(V, v)) {
        V = insert(V, v, size(V)); // add subword j..size
      }
    }
  }
  int nV = size(V);
  intmat GNC[nV][nV]; // graph of n-chains

  // for vertex 1
  for (int i = 2; i <= n + 1; i++) {
    GNC[1,i] = 1; // 1 has an edge to all variables
  }

  // for the other vertices
  for (int i = 2; i <= nV; i++) {
    for (int j = 2; j <= nV; j++) {
      intvec uv = V[i],V[j];

      if (contains(LG, uv)) {
        GNC[i,j] = 1;
      } else {
        // Li p. 177
        // search for a right divisor 'w' of uv in G
        // then check if G doesn't divide the subword uv-1

        // look for a right divisor in LG
        for (int k = 1; k <= size(LG); k++) {
          if (isSF(LG[k], uv)) {
            // w = LG[k]
            if(!ivdivides(LG, uv[1..(size(uv)-1)])) {
              // G doesn't divide uv-1
              GNC[i,j] = 1;
              break;
            }
          }
        }
      }
    }
  }

  return(GNC);
}

static proc contains(list L, def item)
{
  for (int i = 1; i <= size(L); i++) {
    if (L[i] == item) {
      return (1);
    }
  }
  return (0);
}


proc ivDHilbert(list L, int n, list #)
"USAGE: ivDHilbert(L,n[,degbound]); L a list of intmats, n an integer,
@*      degbound an optional integer
RETURN: list
PURPOSE:Computing the K-dimension and the Hilbert series
ASSUME: - basering is a Letterplace ring
@*      - all rows of each intmat correspond to a Letterplace monomial
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds
NOTE: - If L is the list returned, then L[1] is an integer corresponding to the
@*      dimension, L[2] is an intvec which contains the coefficients of the
@*      Hilbert series
@*    - If degbound is set, there will be a degree bound added. By default there
@*      is no degree bound
@*    - n is the number of variables
@*    - If I = L[2] is the intvec returned, then I[k] is the (k-1)-th coefficient of
@*      the Hilbert series.
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example ivDHilbert; shows examples
"
{int degbound = 0;
  if (size(#) > 0){if (typeof(#[1])=="int"){if (#[1] > 0){degbound = #[1];}}}
  checkAssumptions(degbound,L);
  intvec H; int i,dimen;
  H = ivHilbert(L,n,degbound);
  for (i = 1; i <= size(H); i++){dimen = dimen + H[i];}
  L = dimen,H;
  return(L);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  R;
  setring R; // sets basering to Letterplace ring
  //some intmats, which contain monomials in intvec representation as rows
  intmat I1 [2][2] = 1,1,2,2; intmat I2 [1][3]  = 1,2,1;
  intmat J1 [1][2] =  1,1; intmat J2 [2][3] = 2,1,2,1,2,1;
  print(I1);
  print(I2);
  print(J1);
  print(J2);
  list G = I1,I2; // ideal, which is already a Groebner basis
  list I = J1,J2; // ideal, which is already a Groebner basis
  //the procedure without a degree bound
  ivDHilbert(G,2);
  // the procedure with degree bound 5
  ivDHilbert(I,2,5);
}

proc ivDHilbertSickle(list L, int n, list #)
"USAGE: ivDHilbertSickle(L,n[,degbound]); L a list of intmats, n an integer,
@*      degbound an optional integer
RETURN: list
PURPOSE:Computing K-dimension, Hilbert series and mistletoes
ASSUME: - basering is a Letterplace ring.
@*      - All rows of each intmat correspond to a Letterplace monomial.
@*      - If you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If L is the list returned, then L[1] is an integer, L[2] is an intvec
@*      which contains the coefficients of the Hilbert series and L[3]
@*      is a list, containing the mistletoes as intvecs.
@*    - If degbound is set, a degree bound will be added. By default there
@*      is no degree bound.
@*    - n is the number of variables.
@*    - If I = L[2] is the intvec returned, then I[k] is the (k-1)-th
@*      coefficient of the Hilbert series.
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example ivDHilbertSickle; shows examples
"
{int degbound = 0;
  if (size(#) > 0){if (typeof(#[1])=="int"){if (#[1] > 0){degbound = #[1];}}}
  checkAssumptions(degbound,L);
  int i,dimen; list R;
  R = ivSickleHil(L,n,degbound);
  for (i = 1; i <= size(R[1]); i++){dimen = dimen + R[1][i];}
  R[3] = R[2]; R[2] = R[1]; R[1] = dimen;
  return(R);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  R;
  setring R; // sets basering to Letterplace ring
  //some intmats, which contain monomials in intvec representation as rows
  intmat I1 [2][2] = 1,1,2,2; intmat I2 [1][3]  = 1,2,1;
  intmat J1 [1][2] =  1,1; intmat J2 [2][3] = 2,1,2,1,2,1;
  print(I1);
  print(I2);
  print(J1);
  print(J2);
  list G = I1,I2;// ideal, which is already a Groebner basis
  list I = J1,J2; // ideal, which is already a Groebner basis
  ivDHilbertSickle(G,2); // invokes the procedure without a degree bound
  ivDHilbertSickle(I,2,3); // invokes the procedure with degree bound 3
}

proc ivDimCheck(list L, int n)
"USAGE: ivDimCheck(L,n); L a list of intmats, n an integer
RETURN: int, 0 if the dimension is finite, or 1 otherwise
PURPOSE:Decides, whether the K-dimension is finite or not
ASSUME: - basering is a Letterplace ring.
@*      - All rows of each intmat correspond to a Letterplace monomial.
NOTE:   - n is the number of variables.
EXAMPLE: example ivDimCheck; shows examples
"
{checkAssumptions(0,L);
  int i,r;
  intvec P,H;
  for (i = 1; i <= size(L); i++)
  {P[i] = ncols(L[i]);
    if (P[i] == 1) {if (isInMat(H,L[i]) > 0) {ERROR("Quotient algebra is trivial");}}
  }
  if (size(L) == 0) {ERROR("GB is empty, quotient algebra corresponds to free algebra");}
  kill H;
  intmat S; int sd,ld; intvec V;
  sd = P[1]; ld = P[1];
  for (i = 2; i <= size(P); i++)
  {if (P[i] < sd) {sd = P[i];}
    if (P[i] > ld) {ld = P[i];}
  }
  sd = (sd - 1); ld = ld - 1;
  if (ld == 0) { return(allVars(L,P,n));}
  if (sd == 0) { for (i = 1; i <= size(L); i++){if (ncols(L[i]) == 1){S = createStartMat1(n,L[i]); break;}}}
  else {S = createStartMat(sd,n);}
  module M;
  for (i = 1; i <= nrows(S); i++)
  {V = S[i,1..ncols(S)];
    if (findCycle(V,L,P,n,ld,M)) {r = 1; break;}
  }
  return(r);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  R;
  setring R; // sets basering to Letterplace ring
  //some intmats, which contain monomials in intvec representation as rows
  intmat I1 [2][2] = 1,1,2,2; intmat I2 [1][3]  = 1,2,1;
  intmat J1 [1][2] =  1,1; intmat J2 [2][3] = 2,1,2,1,2,1;
  print(I1);
  print(I2);
  print(J1);
  print(J2);
  list G = I1,I2;// ideal, which is already a Groebner basis
  list I = J1,J2; // ideal, which is already a Groebner basis and which
  ivDimCheck(G,2); // invokes the procedure, factor is of finite K-dimension
  ivDimCheck(I,2); // invokes the procedure, factor is not of finite K-dimension
}

proc ivHilbert(list L, int n, list #)
"USAGE: ivHilbert(L,n[,degbound]); L a list of intmats, n an integer,
@*      degbound an optional integer
RETURN: intvec, containing the coefficients of the Hilbert series
PURPOSE:Computing the Hilbert series
ASSUME: - basering is a Letterplace ring.
@*      - all rows of each intmat correspond to a Letterplace monomial
@*      - if you specify a different degree bound degbound,
@*       degbound <= attrib(basering,uptodeg) holds.
NOTE: - If degbound is set, a degree bound  will be added. By default there
@*      is no degree bound.
@*    - n is the number of variables.
@*    - If I is returned, then I[k] is the (k-1)-th coefficient of the Hilbert
@*      series.
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example ivHilbert; shows examples
"
{int degbound = 0;
  if (size(#) > 0) {if (typeof(#[1])=="int"){if (#[1] > 0) {degbound = #[1];}}}
  intvec P,H; int i;
  for (i = 1; i <= size(L); i++)
  {P[i] = ncols(L[i]);
    if (P[i] == 1) {if ( isInMat(H,L[i]) > 0) {ERROR("Quotient algebra is trivial");}}
  }
  if (size(L) == 0) {ERROR("GB is empty, quotient algebra corresponds to free algebra");}
  H[1] = 1;
  checkAssumptions(degbound,L);
  if (degbound == 0)
  {int sd;
    intmat S;
    sd = P[1];
    for (i = 2; i <= size(P); i++) {if (P[i] < sd) {sd = P[i];}}
    sd = (sd - 1);
    if (sd == 0) { for (i = 1; i <= size(L); i++){if (ncols(L[i]) == 1){S = createStartMat1(n,L[i]); break;}}}
    else {S = createStartMat(sd,n);}
    if (intvec(S) == 0) {return(H);}
    for (i = 1; i <= sd; i++) {H = H,(n^i);}
    for (i = 1; i <= nrows(S); i++)
    {intvec St = S[i,1..ncols(S)];
      H = findHCoeff(St,n,L,P,H);
      kill St;
    }
    return(H);
  }
  else
  {for (i = 1; i <= size(P); i++)
    {if (P[i] > degbound) {ERROR("degreebound is too small, GB contains elements of higher degree");}}
    int sd;
    intmat S;
    sd = P[1];
    for (i = 2; i <= size(P); i++) {if (P[i] < sd) {sd = P[i];}}
    sd = (sd - 1);
    if (sd == 0) { for (i = 1; i <= size(L); i++){if (ncols(L[i]) == 1){S = createStartMat1(n,L[i]); break;}}}
    else {S = createStartMat(sd,n);}
    if (intvec(S) == 0) {return(H);}
    for (i = 1; i <= sd; i++) {H = H,(n^i);}
    for (i = 1; i <= nrows(S); i++)
    {intvec St = S[i,1..ncols(S)];
      H = findHCoeff(St,n,L,P,H,degbound);
      kill St;
    }
    return(H);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  R;
  setring R; // sets basering to Letterplace ring
  //some intmats, which contain monomials in intvec representation as rows
  intmat I1 [2][2] = 1,1,2,2; intmat I2 [1][3]  = 1,2,1;
  intmat J1 [1][2] =  1,1; intmat J2 [2][3] = 2,1,2,1,2,1;
  print(I1);
  print(I2);
  print(J1);
  print(J2);
  list G = I1,I2; // ideal, which is already a Groebner basis
  list I = J1,J2; // ideal, which is already a Groebner basis
  ivHilbert(G,2); // invokes the procedure without any degree bound
  ivHilbert(I,2,5); // invokes the procedure with degree bound 5
}


proc ivKDim(list L, int n, list #)
"USAGE: ivKDim(L,n[,degbound]); L a list of intmats,
@*      n an integer, degbound an optional integer
RETURN: int, the K-dimension of A/<L>
PURPOSE:Computing the K-dimension of A/<L>
ASSUME: - basering is a Letterplace ring.
@*      - all rows of each intmat correspond to a Letterplace monomial
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If degbound is set, a degree bound will be added. By default there
@*      is no degree bound.
@*    - n is the number of variables.
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example ivKDim; shows examples
"
{int degbound = 0;
  if (size(#) > 0) {if (typeof(#[1])=="int"){if (#[1] > 0) {degbound = #[1];}}}
  intvec P,H; int i;
  for (i = 1; i <= size(L); i++)
  {P[i] = ncols(L[i]);
    if (P[i] == 1) {if ( isInMat(H,L[i]) > 0) {ERROR("Quotient algebra is trivial");}}
  }
  if (size(L) == 0) {ERROR("GB is empty, quotient algebra corresponds to free algebra");}
  kill H;
  checkAssumptions(degbound,L);
  if (degbound == 0)
  {int sd; int dimen = 1;
    intmat S;
    sd = P[1];
    for (i = 2; i <= size(P); i++) {if (P[i] < sd) {sd = P[i];}}
    sd = (sd - 1);
    if (sd == 0) { for (i = 1; i <= size(L); i++){if (ncols(L[i]) == 1){S = createStartMat1(n,L[i]); break;}}}
    else {S = createStartMat(sd,n);}
    if (intvec(S) == 0) {return(dimen);}
    for (i = 1; i <= sd; i++) {dimen = dimen +(n^i);}
    for (i = 1; i <= nrows(S); i++)
    {intvec St = S[i,1..ncols(S)];
      dimen = dimen + findDimen(St,n,L,P);
      kill St;
    }
    return(dimen);
  }
  else
  {for (i = 1; i <= size(P); i++)
    {if (P[i] > degbound) {ERROR("degreebound is too small, GB contains elements of higher degree");}}
    int sd; int dimen = 1;
    intmat S;
    sd = P[1];
    for (i = 2; i <= size(P); i++) {if (P[i] < sd) {sd = P[i];}}
    sd = (sd - 1);
    if (sd == 0) { for (i = 1; i <= size(L); i++){if (ncols(L[i]) == 1){S = createStartMat1(n,L[i]); break;}}}
    else {S = createStartMat(sd,n);}
    if (intvec(S) == 0) {return(dimen);}
    for (i = 1; i <= sd; i++) {dimen = dimen +(n^i);}
    for (i = 1; i <= nrows(S); i++)
    {intvec St = S[i,1..ncols(S)];
      dimen = dimen + findDimen(St,n,L,P, degbound);
      kill St;
    }
    return(dimen);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  R;
  setring R; // sets basering to Letterplace ring
  //some intmats, which contain monomials in intvec representation as rows
  intmat I1 [2][2] = 1,1,2,2; intmat I2 [1][3]  = 1,2,1;
  intmat J1 [1][2] =  1,1; intmat J2 [2][3] = 2,1,2,1,2,1;
  print(I1);
  print(I2);
  print(J1);
  print(J2);
  list G = I1,I2; // ideal, which is already a Groebner basis
  list I = J1,J2; // ideal, which is already a Groebner basis
  ivKDim(G,2); // invokes the procedure without any degree bound
  ivKDim(I,2,5); // invokes the procedure with degree bound 5
}

proc ivMis2Base(list M)
"USAGE: ivMis2Base(M); M a list of intvecs
RETURN: ideal, a K-base of the given algebra
PURPOSE:Computing the K-base out of given mistletoes
ASSUME: - The mistletoes have to be ordered lexicographically -> OrdMisLex.
@*        Otherwise there might some elements missing.
@*      - basering is a Letterplace ring.
@*      - mistletoes are stored as intvecs, as described in the overview
EXAMPLE: example ivMis2Base; shows examples
"
{
  //checkAssumptions(0,M);
  intvec L,A;
  if (size(M) == 0){ERROR("There are no mistletoes, so it appears your dimension is infinite!");}
  if (isInList(L,M) > 0) {print("1 is a mistletoe, therefore 1 is the only basis element"); return(list(intvec(0)));}
  int i,j,d,s;
  list Rt;
  Rt[1] = intvec(0);
  L = M[1];
  for (i = size(L); 1 <= i; i--) {Rt = insert(Rt,intvec(L[1..i]));}
  for (i = 2; i <= size(M); i++)
  {A = M[i]; L = M[i-1];
    s = size(A);
    if (s > size(L))
    {d = size(L);
      for (j = s; j > d; j--) {Rt = insert(Rt,intvec(A[1..j]));}
      A = A[1..d];
    }
    if (size(L) > s){L = L[1..s];}
    while (A <> L)
    {Rt = insert(Rt, intvec(A));
      if (size(A) > 1)
      {A = A[1..(size(A)-1)];
        L = L[1..(size(L)-1)];
      }
      else {break;}
    }
  }
  return(Rt);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  R;
  setring R; // sets basering to Letterplace ring
  intvec i1 = 1,2; intvec i2 = 2,1,2;
  // the mistletoes are xy and yxy, which are already ordered lexicographically
  list L = i1,i2;
  ivMis2Base(L); // returns the basis of the factor algebra
}


proc ivMis2Dim(list M)
"USAGE: ivMis2Dim(M); M a list of intvecs
RETURN: int, the K-dimension of the given algebra
PURPOSE:Computing the K-dimension out of given mistletoes
ASSUME: - The mistletoes have to be ordered lexicographically -> OrdMisLex.
@*        Otherwise the returned value may differ from the K-dimension.
@*      - basering is a Letterplace ring.
EXAMPLE: example ivMis2Dim; shows examples
"
{checkAssumptions(0,M);
  intvec L;
  if (size(M) == 0){ERROR("There are no mistletoes, so it appears your dimension is infinite!");}
  if (isInList(L,M) > 0) {print("1 is a mistletoe, therefore dim = 1"); return(1);}
  int i,j,d,s;
  j = 1;
  d = 1 + size(M[1]);
  for (i = 1; i < size(M); i++)
  {s = size(M[i]); if (s > size(M[i+1])){s = size(M[i+1]);}
    while ((M[i][j] == M[i+1][j]) && (j <= s)){j = j + 1;}
    d = d + size(M[i+1])- j + 1;
  }
  return(d);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  R;
  setring R; // sets basering to Letterplace ring
  intvec i1 = 1,2; intvec i2 = 2,1,2;
  // the mistletoes are xy and yxy, which are already ordered lexicographically
  list L = i1,i2;
  ivMis2Dim(L); // returns the dimension of the factor algebra
}

proc ivOrdMisLex(list M)
"USAGE: ivOrdMisLex(M); M a list of intvecs
RETURN: list, containing the ordered intvecs of M
PURPOSE:Orders a given set of mistletoes lexicographically
ASSUME: - basering is a Letterplace ring.
- intvecs correspond to monomials
NOTE:   - This is preprocessing, it's not needed if the mistletoes are returned
@*        from the sickle algorithm.
@*      - Each entry of the list returned is an intvec.
EXAMPLE: example ivOrdMisLex; shows examples
"
{checkAssumptions(0,M);
  return(sort(M)[1]);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  intvec i1 = 1,2,1; intvec i2 = 2,2,1; intvec i3 = 1,1; intvec i4 = 2,1,1,1;
  // the corresponding monomials are xyx,y^2x,x^2,yx^3
  list M = i1,i2,i3,i4;
  M;
  ivOrdMisLex(M);// orders the list of monomials
}

proc ivSickle(list L, int n, list #)
"USAGE: ivSickle(L,n,[degbound]); L a list of intmats, n an int, degbound an
@*      optional integer
RETURN: list, containing intvecs, the mistletoes of A/<L>
PURPOSE:Computing the mistletoes for a given Groebner basis L
ASSUME: - basering is a Letterplace ring.
@*      - all rows of each intmat correspond to a Letterplace monomial
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If degbound is set, a degree bound will be added. By default there
@*      is no degree bound.
@*    - n is the number of variables.
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example ivSickle; shows examples
"
{list M;
  int degbound = 0;
  if (size(#) > 0){if (typeof(#[1])=="int"){if (#[1] > 0){degbound = #[1];}}}
  int i;
  intvec P,H;
  for (i = 1; i <= size(L); i++)
  {P[i] = ncols(L[i]);
    if (P[i] == 1) {if (isInMat(H,L[i]) > 0) {ERROR("Quotient algebra is trivial");}}
  }
  if (size(L) == 0) {ERROR("GB is empty, quotient algebra corresponds to free algebra");}
  kill H;
  checkAssumptions(degbound,L);
  if (degbound == 0)
  {intmat S; int sd;
    sd = P[1];
    for (i = 2; i <= size(P); i++) {if (P[i] < sd) {sd = P[i];}}
    sd = (sd - 1);
    if (sd == 0) { for (i = 1; i <= size(L); i++){if (ncols(L[i]) == 1){S = createStartMat1(n,L[i]); break;}}}
    else {S = createStartMat(sd,n);}
    if (intvec(S) == 0) {return(list (intvec(0)));}
    for (i = 1; i <= nrows(S); i++)
    {intvec St = S[i,1..ncols(S)];
      M = M + findmistletoes(St,n,L,P);
      kill St;
    }
    return(M);
  }
  else
  {for (i = 1; i <= size(P); i++)
    {if (P[i] > degbound) {ERROR("degreebound is too small, GB contains elements of higher degree");}}
    intmat S; int sd;
    sd = P[1];
    for (i = 2; i <= size(P); i++) {if (P[i] < sd) {sd = P[i];}}
    sd = (sd - 1);
    if (sd == 0) { for (i = 1; i <= size(L); i++){if (ncols(L[i]) == 1){S = createStartMat1(n,L[i]); break;}}}
    else {S = createStartMat(sd,n);}
    if (intvec(S) == 0) {return(list (intvec(0)));}
    for (i = 1; i <= nrows(S); i++)
    {intvec St = S[i,1..ncols(S)];
      M = M + findmistletoes(St,n,L,P,degbound);
      kill St;
    }
    return(M);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  //some intmats, which contain monomials in intvec representation as rows
  intmat I1 [2][2] = 1,1,2,2; intmat I2 [1][3]  = 1,2,1;
  intmat J1 [1][2] =  1,1; intmat J2 [2][3] = 2,1,2,1,2,1;
  print(I1);
  print(I2);
  print(J1);
  print(J2);
  list G = I1,I2; // ideal, which is already a Groebner basis
  list I =  J1,J2; // ideal, which is already a Groebner basis
  ivSickle(G,2); // invokes the procedure without any degree bound
  ivSickle(I,2,5); // invokes the procedure with degree bound 5
}

proc ivSickleDim(list L, int n, list #)
"USAGE: ivSickleDim(L,n[,degbound]); L a list of intmats, n an integer, degbound
@*      an optional integer
RETURN: list
PURPOSE:Computing mistletoes and the K-dimension
ASSUME: - basering is a Letterplace ring.
@*      - all rows of each intmat correspond to a Letterplace monomial
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If L is the list returned, then L[1] is an integer, L[2] is a list,
@*      containing the mistletoes as intvecs.
@*    - If degbound is set, a degree bound will be added. By default there
@*      is no degree bound.
@*    - n is the number of variables.
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example ivSickleDim; shows examples
"
{list M;
  int degbound = 0;
  if (size(#) > 0){if (typeof(#[1])=="int"){if (#[1] > 0){degbound = #[1];}}}
  int i,dimen; list R;
  intvec P,H;
  for (i = 1; i <= size(L); i++)
  {P[i] = ncols(L[i]);
    if (P[i] == 1) {if (isInMat(H,L[i]) > 0) {ERROR("Quotient algebra is trivial, dimension equals zero");}}
  }
  if (size(L) == 0) {ERROR("GB is empty, quotient algebra corresponds to free algebra");}
  kill H;
  checkAssumptions(degbound,L);
  if (degbound == 0)
  {int sd; dimen = 1;
    intmat S;
    sd = P[1];
    for (i = 2; i <= size(P); i++) {if (P[i] < sd) {sd = P[i];}}
    sd = (sd - 1);
    if (sd == 0) { for (i = 1; i <= size(L); i++){if (ncols(L[i]) == 1){S = createStartMat1(n,L[i]); break;}}}
    else {S = createStartMat(sd,n);}
    if (intvec(S) == 0) {return(list(dimen,list(intvec(0))));}
    for (i = 1; i <= sd; i++) {dimen = dimen +(n^i);}
    R[1] = dimen;
    for (i = 1; i <= nrows(S); i++)
    {intvec St = S[i,1..ncols(S)];
      R = findMisDim(St,n,L,P,R);
      kill St;
    }
    return(R);
  }
  else
  {for (i = 1; i <= size(P); i++)
    {if (P[i] > degbound) {ERROR("degreebound is too small, GB contains elements of higher degree");}}
    int sd; dimen = 1;
    intmat S;
    sd = P[1];
    for (i = 2; i <= size(P); i++) {if (P[i] < sd) {sd = P[i];}}
    sd = (sd - 1);
    if (sd == 0) { for (i = 1; i <= size(L); i++){if (ncols(L[i]) == 1){S = createStartMat1(n,L[i]); break;}}}
    else {S = createStartMat(sd,n);}
    if (intvec(S) == 0) {return(list(dimen,list(intvec(0))));}
    for (i = 1; i <= sd; i++) {dimen = dimen +(n^i);}
    R[1] = dimen;
    for (i = 1; i <= nrows(S); i++)
    {intvec St = S[i,1..ncols(S)];
      R = findMisDim(St,n,L,P,R,degbound);
      kill St;
    }
    return(R);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  //some intmats, which contain monomials in intvec representation as rows
  intmat I1 [2][2] = 1,1,2,2; intmat I2 [1][3]  = 1,2,1;
  intmat J1 [1][2] =  1,1; intmat J2 [2][3] = 2,1,2,1,2,1;
  print(I1);
  print(I2);
  print(J1);
  print(J2);
  list G = I1,I2;// ideal, which is already a Groebner basis
  list I =  J1,J2; // ideal, which is already a Groebner basis
  ivSickleDim(G,2); // invokes the procedure without any degree bound
  ivSickleDim(I,2,5); // invokes the procedure with degree bound 5
}

proc ivSickleHil(list L, int n, list #)
"USAGE:ivSickleHil(L,n[,degbound]); L a list of intmats, n an integer,
@*     degbound an optional integer
RETURN: list
PURPOSE:Computing the mistletoes and the Hilbert series
ASSUME: - basering is a Letterplace ring.
@*      - all rows of each intmat correspond to a Letterplace monomial
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If L is the list returned, then L[1] is an intvec, L[2] is a list,
@*      containing the mistletoes as intvecs.
@*    - If degbound is set, a degree bound will be added. By default there
@*      is no degree bound.
@*    - n is the number of variables.
@*    - If I = L[1] is the intvec returned, then I[k] is the (k-1)-th
@*      coefficient of the Hilbert series.
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example ivSickleHil; shows examples
"
{int degbound = 0;
  if (size(#) > 0) {if (typeof(#[1])=="int"){if (#[1] > 0) {degbound = #[1];}}}
  intvec P,H; int i; list R;
  for (i = 1; i <= size(L); i++)
  {P[i] = ncols(L[i]);
    if (P[i] == 1) {if ( isInMat(H,L[i]) > 0) {ERROR("Quotient algebra is trivial");}}
  }
  if (size(L) == 0) {ERROR("GB is empty, quotient algebra corresponds to free algebra");}
  H[1] = 1;
  checkAssumptions(degbound,L);
  if (degbound == 0)
  {int sd;
    intmat S;
    sd = P[1];
    for (i = 2; i <= size(P); i++) {if (P[i] < sd) {sd = P[i];}}
    sd = (sd - 1);
    if (sd == 0) { for (i = 1; i <= size(L); i++){if (ncols(L[i]) == 1){S = createStartMat1(n,L[i]); break;}}}
    else {S = createStartMat(sd,n);}
    if (intvec(S) == 0) {return(list(H,list(intvec (0))));}
    for (i = 1; i <= sd; i++) {H = H,(n^i);}
    R[1] = H; kill H;
    for (i = 1; i <= nrows(S); i++)
    {intvec St = S[i,1..ncols(S)];
      R = findHCoeffMis(St,n,L,P,R);
      kill St;
    }
    return(R);
  }
  else
  {for (i = 1; i <= size(P); i++)
    {if (P[i] > degbound) {ERROR("degreebound is too small, GB contains elements of higher degree");}}
    int sd;
    intmat S;
    sd = P[1];
    for (i = 2; i <= size(P); i++) {if (P[i] < sd) {sd = P[i];}}
    sd = (sd - 1);
    if (sd == 0) { for (i = 1; i <= size(L); i++){if (ncols(L[i]) == 1){S = createStartMat1(n,L[i]); break;}}}
    else {S = createStartMat(sd,n);}
    if (intvec(S) == 0) {return(list(H,list(intvec(0))));}
    for (i = 1; i <= sd; i++) {H = H,(n^i);}
    R[1] = H; kill H;
    for (i = 1; i <= nrows(S); i++)
    {intvec St = S[i,1..ncols(S)];
      R = findHCoeffMis(St,n,L,P,R,degbound);
      kill St;
    }
    return(R);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  //some intmats, which contain monomials in intvec representation as rows
  intmat I1[2][2] = 1,1,2,2; intmat I2[1][3]  = 1,2,1;
  intmat J1[1][2] =  1,1; intmat J2[2][3] = 2,1,2,1,2,1;
  print(I1);
  print(I2);
  print(J1);
  print(J2);
  list G = I1,I2;// ideal, which is already a Groebner basis
  list I =  J1,J2; // ideal, which is already a Groebner basis
  ivSickleHil(G,2); // invokes the procedure without any degree bound
  ivSickleHil(I,2,5); // invokes the procedure with degree bound 5
}

proc lpDHilbert(ideal G, list #)
"USAGE: lpDHilbert(G[,degbound,n]); G an ideal, degbound, n optional integers
RETURN: list
PURPOSE:Computing K-dimension and Hilbert series, starting with a lp-ideal
ASSUME: - basering is a Letterplace ring.
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If L is the list returned, then L[1] is an integer corresponding to the
@*      dimension, L[2] is an intvec which contains the coefficients of the
@*      Hilbert series
@*    - If degbound is set, there will be a degree bound added. 0 means no
@*      degree bound. Default: attrib(basering,uptodeg).
@*    - n can be set to a different number of variables.
@*      Default: n = attrib(basering, lV).
@*    - If I = L[2] is the intvec returned, then I[k] is the (k-1)-th
@*      coefficient of the Hilbert series.
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example lpDHilbert; shows examples
"
{int degbound = attrib(basering,"uptodeg");int n = attrib(basering, "lV");
  if (size(#) > 0){if (typeof(#[1])=="int"){if (#[1] >= 0){degbound = #[1];}}}
  if (size(#) > 1){if (typeof(#[1])=="int"){if (#[2] > 0){n = #[2];}}}
  list L;
  L = lp2ivId(normalize(lead(G)));
  return(ivDHilbert(L,n,degbound));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal G = x(1)*x(2), y(1)*y(2),x(1)*y(2)*x(3); // ideal G contains a
  //Groebner basis
  lpDHilbert(G,5,2); // invokes procedure with degree bound 5 and 2 variables
  // note that the optional parameters are not necessary, due to the finiteness
  // of the K-dimension of the factor algebra
  lpDHilbert(G); // procedure with ring parameters
  lpDHilbert(G,0); // procedure without degreebound
}

proc lpDHilbertSickle(ideal G, list #)
"USAGE: lpDHilbertSickle(G[,degbound,n]); G an ideal, degbound, n optional
@*      integers
RETURN: list
PURPOSE:Computing K-dimension, Hilbert series and mistletoes at once
ASSUME: - basering is a Letterplace ring.
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If L is the list returned, then L[1] is an integer, the K-dimension,
@*      L[2] is an intvec, the Hilbert series and L[3] is an ideal,
@*      the mistletoes
@*    - If degbound is set, there will be a degree bound added. 0 means no
@*      degree bound. Default: attrib(basering,uptodeg).
@*    - n can be set to a different number of variables.
@*      Default: n = attrib(basering, lV).
@*    - If I = L[1] is the intvec returned, then I[k] is the (k-1)-th
@*      coefficient of the Hilbert series.
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example lpDHilbertSickle; shows examples
"
{int degbound = attrib(basering,"uptodeg");int n = attrib(basering, "lV");
  if (size(#) > 0){if (typeof(#[1])=="int"){if (#[1] >= 0){degbound = #[1];}}}
  if (size(#) > 1){if (typeof(#[1])=="int"){if (#[2] > 0){n = #[2];}}}
  list L;
  L = lp2ivId(normalize(lead(G)));
  L = ivDHilbertSickle(L,n,degbound);
  L[3] =  ivL2lpI(L[3]);
  return(L);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal G = x(1)*x(2), y(1)*y(2),x(1)*y(2)*x(3); // ideal G contains a
  //Groebner basis
  lpDHilbertSickle(G,5,2); //invokes procedure with degree bound 5 and 2 variables
  // note that the optional parameters are not necessary, due to the finiteness
  // of the K-dimension of the factor algebra
  lpDHilbertSickle(G); // procedure with ring parameters
  lpDHilbertSickle(G,0); // procedure without degreebound
}

proc lpHilbert(ideal G, list #)
"USAGE: lpHilbert(G[,degbound,n]); G an ideal, degbound, n optional integers
RETURN: intvec, containing the coefficients of the Hilbert series
PURPOSE:Computing the Hilbert series
ASSUME: - basering is a Letterplace ring.
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If degbound is set, there will be a degree bound added. 0 means no
@*      degree bound. Default: attrib(basering,uptodeg).
@*    - n is the number of variables, which can be set to a different number.
@*      Default: attrib(basering, lV).
@*    - If I is returned, then I[k] is the (k-1)-th coefficient of the Hilbert
@*      series.
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example lpHilbert; shows examples
"
{int degbound = attrib(basering,"uptodeg");int n = attrib(basering, "lV");
  if (size(#) > 0){if (typeof(#[1])=="int"){if (#[1] >= 0){degbound = #[1];}}}
  if (size(#) > 1){if (typeof(#[1])=="int"){if (#[2] > 0){n = #[2];}}}
  list L;
  L = lp2ivId(normalize(lead(G)));
  return(ivHilbert(L,n,degbound));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal G = x(1)*x(2), y(1)*y(2),x(1)*y(2)*x(3); // ideal G contains a
  //Groebner basis
  lpHilbert(G,5,2); // invokes procedure with degree bound 5 and 2 variables
  // note that the optional parameters are not necessary, due to the finiteness
  // of the K-dimension of the factor algebra
  lpDHilbert(G); // procedure with ring parameters
  lpDHilbert(G,0); // procedure without degreebound
}

proc lpDimCheck(ideal G)
"USAGE: lpDimCheck(G);
RETURN: int, 1 if K-dimension of the factor algebra is infinite, 0 otherwise
PURPOSE:Checking a factor algebra for finiteness of the K-dimension
ASSUME: - basering is a Letterplace ring.
EXAMPLE: example lpDimCheck; shows examples
"
{int n = attrib(basering,"lV");
  list L;
  ideal R;
  R = normalize(lead(G));
  L = lp2ivId(R);
  return(ivDimCheck(L,n));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal G = x(1)*x(2), y(1)*y(2),x(1)*y(2)*x(3);
  // Groebner basis
  ideal I = x(1)*x(2), y(1)*x(2)*y(3), x(1)*y(2)*x(3);
  // Groebner basis
  lpDimCheck(G); // invokes procedure, factor algebra is of finite K-dimension
  lpDimCheck(I); // invokes procedure, factor algebra is of infinite Kdimension
}

proc lpKDim(ideal G, list #)
"USAGE: lpKDim(G[,degbound, n]); G an ideal, degbound, n optional integers
RETURN: int, the K-dimension of the factor algebra
PURPOSE:Computing the K-dimension of a factor algebra, given via an ideal
ASSUME: - basering is a Letterplace ring
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If degbound is set, there will be a degree bound added. 0 means no
@*      degree bound. Default: attrib(basering, uptodeg).
@*    - n is the number of variables, which can be set to a different number.
@*      Default: attrib(basering, lV).
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example lpKDim; shows examples
"
{int degbound = attrib(basering, "uptodeg");int n = attrib(basering, "lV");
  if (size(#) > 0){if (typeof(#[1])=="int"){if (#[1] >= 0){degbound = #[1];}}}
  if (size(#) > 1){if (typeof(#[1])=="int"){if (#[2] > 0){n = #[2];}}}
  list L;
  L = lp2ivId(normalize(lead(G)));
  return(ivKDim(L,n,degbound));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal G = x(1)*x(2), y(1)*y(2),x(1)*y(2)*x(3);
  // ideal G contains a Groebner basis
  lpKDim(G); //procedure invoked with ring parameters
  // the factor algebra is finite, so the degree bound given by the Letterplace
  // ring is not necessary
  lpKDim(G,0); // procedure without any degree bound
}

proc lpMis2Base(ideal M)
"USAGE: lpMis2Base(M); M an ideal
RETURN: ideal, a K-basis of the factor algebra
PURPOSE:Computing a K-basis out of given mistletoes
ASSUME: - basering is a Letterplace ring. G is a Letterplace ideal.
@*      - M contains only monomials
NOTE:   - The mistletoes have to be ordered lexicographically -> OrdMisLex.
EXAMPLE: example lpMis2Base; shows examples
"
{list L;
  L = lpId2ivLi(M);
  return(ivL2lpI(ivMis2Base(L)));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal L = x(1)*y(2),y(1)*x(2)*y(3);
  // ideal containing the mistletoes
  lpMis2Base(L); // returns the K-basis of the factor algebra
}

proc lpMis2Dim(ideal M)
"USAGE: lpMis2Dim(M); M an ideal
RETURN: int, the K-dimension of the factor algebra
PURPOSE:Computing the K-dimension out of given mistletoes
ASSUME: - basering is a Letterplace ring.
@*      - M contains only monomials
NOTE:   - The mistletoes have to be ordered lexicographically -> OrdMisLex.
EXAMPLE: example lpMis2Dim; shows examples
"
{list L;
  L = lpId2ivLi(M);
  return(ivMis2Dim(L));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal L = x(1)*y(2),y(1)*x(2)*y(3);
  // ideal containing the mistletoes
  lpMis2Dim(L); // returns the K-dimension of the factor algebra
}

proc lpOrdMisLex(ideal M)
"USAGE: lpOrdMisLex(M); M an ideal of mistletoes
RETURN: ideal, containing the mistletoes, ordered lexicographically
PURPOSE:A given set of mistletoes is ordered lexicographically
ASSUME: - basering is a Letterplace ring.
NOTE:   This is preprocessing, it is not needed if the mistletoes are returned
@*      from the sickle algorithm.
EXAMPLE: example lpOrdMisLex; shows examples
"
{return(ivL2lpI(sort(lpId2ivLi(M))[1]));}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal M = x(1)*y(2)*x(3), y(1)*y(2)*x(3), x(1)*x(2), y(1)*x(2)*x(3)*x(4);
  // some monomials
  lpOrdMisLex(M); // orders the monomials lexicographically
}

proc lpSickle(ideal G,  list #)
"USAGE: lpSickle(G[,degbound,n]); G an ideal, degbound, n optional integers
RETURN: ideal
PURPOSE:Computing the mistletoes of K[X]/<G>
ASSUME: - basering is a Letterplace ring.
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If degbound is set, there will be a degree bound added. 0 means no
@*      degree bound. Default: attrib(basering,uptodeg).
@*    - n is the number of variables, which can be set to a different number.
@*      Default: attrib(basering, lV).
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example lpSickle; shows examples
"
{int degbound = attrib(basering,"uptodeg"); int n = attrib(basering, "lV");
  if (size(#) > 0){if (typeof(#[1])=="int"){if (#[1] >= 0){degbound = #[1];}}}
  if (size(#) > 1){if (typeof(#[1])=="int"){if (#[2] > 0){n = #[2];}}}
  list L; ideal R;
  R = normalize(lead(G));
  L = lp2ivId(R);
  L = ivSickle(L,n,degbound);
  R = ivL2lpI(L);
  return(R);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal G = x(1)*x(2), y(1)*y(2),x(1)*y(2)*x(3); // ideal G contains a
  //Groebner basis
  lpSickle(G); //invokes the procedure with ring parameters
  // the factor algebra is finite, so the degree bound given by the Letterplace
  // ring is not necessary
  lpSickle(G,0); // procedure without any degree bound
}

proc lpSickleDim(ideal G, list #)
"USAGE: lpSickleDim(G[,degbound,n]); G an ideal, degbound, n optional integers
RETURN: list
PURPOSE:Computing the K-dimension and the mistletoes
ASSUME: - basering is a Letterplace ring.
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If L is the list returned, then L[1] is an integer, the K-dimension,
@*      L[2] is an ideal, the mistletoes.
@*    - If degbound is set, there will be a degree bound added. 0 means no
@*      degree bound. Default: attrib(basering,uptodeg).
@*    - n is the number of variables, which can be set to a different number.
@*      Default: attrib(basering, lV).
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example lpSickleDim; shows examples
"
{int degbound = attrib(basering,"uptodeg");int n = attrib(basering, "lV");
  if (size(#) > 0){if (typeof(#[1])=="int"){if (#[1] >= 0){degbound = #[1];}}}
  if (size(#) > 1){if (typeof(#[1])=="int"){if (#[2] > 0){n = #[2];}}}
  list L;
  L = lp2ivId(normalize(lead(G)));
  L = ivSickleDim(L,n,degbound);
  L[2] = ivL2lpI(L[2]);
  return(L);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal G = x(1)*x(2), y(1)*y(2),x(1)*y(2)*x(3); // ideal G contains a
  //Groebner basis
  lpSickleDim(G); // invokes the procedure with ring parameters
  // the factor algebra is finite, so the degree bound given by the Letterplace
  // ring is not necessary
  lpSickleDim(G,0); // procedure without any degree bound
}

proc lpSickleHil(ideal G, list #)
"USAGE: lpSickleHil(G);
RETURN: list
PURPOSE:Computing the Hilbert series and the mistletoes
ASSUME: - basering is a Letterplace ring.
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE: - If L is the list returned, then L[1] is an intvec, corresponding to the
@*      Hilbert series, L[2] is an ideal, the mistletoes.
@*    - If degbound is set, there will be a degree bound added. 0 means no
@*      degree bound. Default: attrib(basering,uptodeg).
@*    - n is the number of variables, which can be set to a different number.
@*      Default: attrib(basering, lV).
@*    - If I = L[1] is the intvec returned, then I[k] is the (k-1)-th
@*      coefficient of the Hilbert series.
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example lpSickleHil; shows examples
"
{int degbound = attrib(basering,"uptodeg");int n = attrib(basering, "lV");
  if (size(#) > 0){if (typeof(#[1])=="int"){if (#[1] >= 0){degbound = #[1];}}}
  if (size(#) > 1){if (typeof(#[1])=="int"){if (#[2] > 0){n = #[2];}}}
  list L;
  L = lp2ivId(normalize(lead(G)));
  L = ivSickleHil(L,n,degbound);
  L[2] =  ivL2lpI(L[2]);
  return(L);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal G = x(1)*x(2), y(1)*y(2),x(1)*y(2)*x(3); // ideal G contains a
  //Groebner basis
  lpSickleHil(G); // invokes the procedure with ring parameters
  // the factor algebra is finite, so the degree bound given by the Letterplace
  // ring is not necessary
  lpSickleHil(G,0); // procedure without any degree bound
}

proc sickle(ideal G, list #)
"USAGE: sickle(G[,m, d, h, degbound]); G an ideal; m,d,h,degbound optional
@*      integers
RETURN: list
PURPOSE:Allowing the user to access all procs with one command
ASSUME: - basering is a Letterplace ring.
@*      - if you specify a different degree bound degbound,
@*        degbound <= attrib(basering,uptodeg) holds.
NOTE:   The returned object will always be a list, but the entries of the
@*      returned list may be very different
@* case m=1,d=1,h=1: see lpDHilbertSickle
@* case m=1,d=1,h=0: see lpSickleDim
@* case m=1,d=0,h=1: see lpSickleHil
@* case m=1,d=0,h=0: see lpSickle (this is the default case)
@* case m=0,d=1,h=1: see lpDHilbert
@* case m=0,d=1,h=0: see lpKDim
@* case m=0,d=0,h=1: see lpHilbert
@* case m=0,d=0,h=0: returns an error
@*    - If degbound is set, there will be a degree bound added. 0 means no
@*      degree bound. Default: attrib(basering,uptodeg).
@*    - If the K-dimension is known to be infinite, a degree bound is needed
EXAMPLE: example sickle; shows examples
"
{int m,d,h,degbound;
  m = 1; d = 0; h = 0; degbound = attrib(basering,"uptodeg");
  if (size(#) > 0) {if (typeof(#[1])=="int"){if (#[1] < 1) {m = 0;}}}
  if (size(#) > 1) {if (typeof(#[1])=="int"){if (#[2] > 0) {d = 1;}}}
  if (size(#) > 2) {if (typeof(#[1])=="int"){if (#[3] > 0) {h = 1;}}}
  if (size(#) > 3) {if (typeof(#[1])=="int"){if (#[4] >= 0) {degbound = #[4];}}}
  if (m == 1)
  {if (d == 0)
    {if (h == 0) {return(lpSickle(G,degbound,attrib(basering,"lV")));}
      else        {return(lpSickleHil(G,degbound,attrib(basering,"lV")));}
    }
    else
    {if (h == 0) {return(lpSickleDim(G,degbound,attrib(basering,"lV")));}
      else {return(lpDHilbertSickle(G,degbound,attrib(basering,"lV")));}
    }
  }
  else
  {if (d == 0)
    {if (h == 0) {ERROR("You request to do nothing, so relax and do so");}
      else        {return(lpHilbert(G,degbound,attrib(basering,"lV")));}
    }
    else
    {if (h == 0) {return(lpKDim(G,degbound,attrib(basering,"lV")));}
      else {return(lpDHilbert(G,degbound,attrib(basering,"lV")));}
    }
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal G = x(1)*x(2), y(1)*y(2),x(1)*y(2)*x(3);
  // G contains a Groebner basis
  sickle(G,1,1,1); // computes mistletoes, K-dimension and the Hilbert series
  sickle(G,1,0,0); // computes mistletoes only
  sickle(G,0,1,0); // computes K-dimension only
  sickle(G,0,0,1); // computes Hilbert series only
}

///////////////////////////////////////////////////////////////////////////////
/* vl: new stuff for conversion to Magma and to SD
todo: doc, example
 */
proc extractVars(r)
{
  int i = 1;
  int j = 1;
  string candidate;
  list result = list();
  for (i = 1; i<=nvars(r);i++)
  {
    candidate = string(var(i))[1,find(string(var(i)),"(")-1];
    if (!inList(result, candidate))
    {
      result = insert(result,candidate,size(result));
    }
  }
  return(result);
}

proc letterPlacePoly2MagmaString(poly h)
{
  int pos;
  string s = string(h);
  while(find(s,"("))
  {
    pos = find(s,"(");
    while(s[pos]!=")")
    {
      s = s[1,pos-1]+s[pos+1,size(s)-pos];
    }
    if (size(s)!=pos)
    {
      s = s[1,pos-1]+s[pos+1,size(s)-pos]; // The last (")")
    }
    else
    {
      s = s[1,pos-1];
    }
  }
  return(s);
}

proc letterPlaceIdeal2SD(ideal I, int upToDeg)
{
  int i;
  print("Don't forget to fill in the formal Data in the file");
  string result = "<?xml version=\"1.0\"?>"+newline+"<FREEALGEBRA createdAt=\"\" createdBy=\"Singular\" id=\"FREEALGEBRA/\">"+newline;
  result = result + "<vars>"+string(extractVars(basering))+"</vars>"+newline;
  result = result + "<basis>"+newline;
  for (i = 1;i<=size(I);i++)
  {
    result = result + "<poly>"+letterPlacePoly2MagmaString(I[i])+"</poly>"+newline;
  }
  result = result + "</basis>"+newline;
  result = result + "<uptoDeg>"+ string(upToDeg)+"</uptoDeg>"+newline;
  result = result + "<Comment></Comment>"+newline;
  result = result + "<Version></Version>"+newline;
  result = result + "</FREEALGEBRA>";
  return(result);
}


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


proc tst_fpadim()
{
  example ivDHilbert;
  example ivDHilbertSickle;
  example ivDimCheck;
  example ivHilbert;
  example ivKDim;
  example ivMis2Base;
  example ivMis2Dim;
  example ivOrdMisLex;
  example ivSickle;
  example ivSickleHil;
  example ivSickleDim;
  example lpDHilbert;
  example lpDHilbertSickle;
  example lpHilbert;
  example lpDimCheck;
  example lpKDim;
  example lpMis2Base;
  example lpMis2Dim;
  example lpOrdMisLex;
  example lpSickle;
  example lpSickleHil;
  example lpSickleDim;
  example sickle;
  example ivL2lpI;
  example iv2lp;
  example iv2lpList;
  example iv2lpMat;
  example lp2iv;
  example lp2ivId;
  example lpId2ivLi;
  example lpGkDim;
  example lpGlDimBound;
  example lpSubstitute;
}


/*proc lpSubstituteExpandRing(poly f, list s1, list s2) {*/
/*int minDegBound = calcSubstDegBound(f,s1,s2);*/
/**/
/*def R = basering; // curr lp ring*/
/*setring ORIGINALRING; // non lp ring TODO*/
/*def R1 = makeLetterplaceRing(minDegBound);*/
/*setring R1;*/
/**/
/*poly g = lpSubstitute(imap(R,f), imap(R,s1), imap(R,s2));*/
/**/
/*return (R1); // return the new ring*/
/*}*/

proc lpSubstitute(poly f, ideal s1, ideal s2, list #)
"USAGE: lpSubstitute(f,s1,s2[,G]); f letterplace polynomial, s1 list (ideal) of variables
@*      to replace, s2 list (ideal) of polynomials to replace with, G optional ideal to
@*      reduce with.
RETURN: poly, the substituted polynomial
ASSUME: - basering is a Letterplace ring
@*      - G is a groebner basis,
@*      - the current ring has a sufficient degbound (can be calculated with
    @*    calcSubstDegBound())
EXAMPLE: example lpSubstitute; shows examples
"
{
  ideal G;
  if (size(#) > 0) {
    if (typeof(#[1])=="ideal") {
      G = #[1];
    }
  }

  poly fs;
  for (int i = 1; i <= size(f); i++) {
    poly fis = leadcoef(f[i]);
    intvec ivfi = lp2iv(f[i]);
    for (int j = 1; j <= size(ivfi); j++) {
      int varindex = ivfi[j];
      int subindex = lpIndexOf(s1, var(varindex));
      if (subindex > 0) {
        s2[subindex] = lpNF(s2[subindex],G);
        fis = lpMult(fis, s2[subindex]);
      } else {
        fis = lpMult(fis, lpNF(iv2lp(varindex),G));
      }
      /*fis = lpNF(fis,G);*/
    }
    fs = fs + fis;
  }
  fs = lpNF(fs, G);
  return (fs);
}
example {
  LIB "fpadim.lib";
  ring r = 0,(x,y,z),dp;
  def R = makeLetterplaceRing(4);
  setring R;

  ideal G = x(1)*y(2); // optional

  poly f = 3*x(1)*x(2)+y(1)*x(2);
  ideal s1 = x(1), y(1);
  ideal s2 = y(1)*z(2)*z(3), x(1);

  // the substitution probably needs a higher degbound
  int minDegBound = calcSubstDegBounds(f,s1,s2);
  setring r;
  def R1 = makeLetterplaceRing(minDegBound);
  setring R1;

  // the last parameter is optional
  lpSubstitute(imap(R,f), imap(R,s1), imap(R,s2), imap(R,G));
}
example {
  LIB "fpadim.lib";
  ring r = 0,(x,y,z),dp;
  def R = makeLetterplaceRing(4);
  setring R;

  poly f = 3*x(1)*x(2)+y(1)*x(2);
  poly g = z(1)*x(2)+y(1);
  poly h = 7*x(1)*z(2)+x(1);
  ideal I = f,g,h;
  ideal s1 = x(1), y(1);
  ideal s2 = y(1)*z(2)*z(3), x(1);

  int minDegBound = calcSubstDegBounds(I,s1,s2);
  setring r;
  def R1 = makeLetterplaceRing(minDegBound);
  setring R1;

  ideal I = imap(R,I);
  ideal s1 = imap(R,s1);
  ideal s2 = imap(R,s2);
  for (int i = 1; i <= size(I); i++) {
    lpSubstitute(I[i], s1, s2);
  }
}

static proc lpIndexOf(ideal I, poly p) {
  for (int i = 1; i <= size(I); i++) {
    if (I[i] == p) {
      return (i);
    }
  }
  return (-1);
}

static proc ivIndexOf(list L, intvec iv) {
  for (int i = 1; i <= size(L); i++) {
    if (L[i] == iv) {
      return (i);
    }
  }
  return (-1);
}


proc calcSubstDegBound(poly f, ideal s1, ideal s2)
"USAGE: calcSubstDegBound(f,s1,s2); f letterplace polynomial, s1 list (ideal) of variables
@*      to replace, s2 list (ideal) of polynomials to replace with
RETURN: int, the min degbound required to perform the substitution
ASSUME: - basering is a Letterplace ring
EXAMPLE: example lpSubstitute; shows examples
"
{
  int maxDegBound = 0;
  for (int i = 1; i <= size(f); i++) {
    intvec ivfi = lp2iv(f[i]);
    int tmpDegBound;
    for (int j = 1; j <= size(ivfi); j++) {
      int varindex = ivfi[j];
      int subindex = lpIndexOf(s1, var(varindex));
      if (subindex > 0) {
        tmpDegBound = tmpDegBound + deg(s2[subindex]);
      } else {
        tmpDegBound = tmpDegBound + 1;
      }
    }
    if (tmpDegBound > maxDegBound) {
      maxDegBound = tmpDegBound;
    }
  }

  // increase degbound by 50% when ideal is provided
  // needed for lpNF
  maxDegBound = maxDegBound + maxDegBound/2;

  return (maxDegBound);
}

// convenience method
proc calcSubstDegBounds(ideal I, ideal s1, ideal s2)
"USAGE: calcSubstDegBounds(I,s1,s2); I list (ideal) of letterplace polynomials, s1 list (ideal)
@*      of variables to replace, s2 list (ideal) of polynomials to replace with
RETURN: int, the min degbound required to perform all of the substitutions
ASSUME: - basering is a Letterplace ring
EXAMPLE: example lpSubstitute; shows examples
"
{
  int maxDegBound = 0;
  for (int i = 1; i <= size(I); i++) {
    int tmpDegBound = calcSubstDegBound(I[i], s1, s2, #);
    if (tmpDegBound > maxDegBound) {
      maxDegBound = tmpDegBound;
    }
  }
  return (maxDegBound);
}


/*
   Here are some examples one may try. Just copy them into your console.
   These are relations for braid groups, up to degree d:


   LIB "fpadim.lib";
   ring r = 0,(x,y,z),dp;
   int d =10; // degree
   def R = makeLetterplaceRing(d);
   setring R;
   ideal I = y(1)*x(2)*y(3) - z(1)*y(2)*z(3), x(1)*y(2)*x(3) - z(1)*x(2)*y(3),
   z(1)*x(2)*z(3) - y(1)*z(2)*x(3), x(1)*x(2)*x(3) + y(1)*y(2)*y(3) +
   z(1)*z(2)*z(3) + x(1)*y(2)*z(3);
   option(prot);
   option(redSB);option(redTail);option(mem);
   ideal J = system("freegb",I,d,3);
   lpDimCheck(J);
   sickle(J,1,1,1,d);//Computes mistletoes, K-dimension and the Hilbert series



   LIB "fpadim.lib";
   ring r = 0,(x,y,z),dp;
   int d =11; // degree
   def R = makeLetterplaceRing(d);
   setring R;
   ideal I = y(1)*x(2)*y(3) - z(1)*y(2)*z(3), x(1)*y(2)*z(3) - z(1)*x(2)*y(3),
   z(1)*x(2)*z(3) - y(1)*z(2)*x(3), x(1)*x(2)*x(3) + y(1)*y(2)*y(3) +
   z(1)*z(2)*z(3) + x(1)*y(2)*z(3);
   option(prot);
   option(redSB);option(redTail);option(mem);
   ideal J = system("freegb",I,d,3);
   lpDimCheck(J);
   sickle(J,1,1,1,d);



   LIB "fpadim.lib";
   ring r = 0,(x,y,z),dp;
   int d  = 6; // degree
   def R  = makeLetterplaceRing(d);
   setring R;
   ideal I = y(1)*x(2)*y(3) - z(1)*y(2)*z(3), x(1)*y(2)*x(3) - z(1)*x(2)*y(3),
   z(1)*x(2)*z(3) - y(1)*z(2)*x(3), x(1)*x(2)*x(3) -2*y(1)*y(2)*y(3) + 3*z(1)*z(2)*z(3) -4*x(1)*y(2)*z(3) + 5*x(1)*z(2)*z(3)- 6*x(1)*y(2)*y(3) +7*x(1)*x(2)*z(3) - 8*x(1)*x(2)*y(3);
   option(prot);
   option(redSB);option(redTail);option(mem);
   ideal J = system("freegb",I,d,3);
   lpDimCheck(J);
   sickle(J,1,1,1,d);
 */

/*
   Here are some examples, which can also be found in [studzins]:

// takes up to 880Mb of memory
LIB "fpadim.lib";
ring r = 0,(x,y,z),dp;
int d =10; // degree
def R = makeLetterplaceRing(d);
setring R;
ideal I =
z(1)*z(2)*z(3)*z(4) + y(1)*x(2)*y(3)*x(4) - x(1)*y(2)*y(3)*x(4) - 3*z(1)*y(2)*x(3)*z(4), x(1)*x(2)*x(3) + y(1)*x(2)*y(3) - x(1)*y(2)*x(3), z(1)*y(2)*x(3)-x(1)*y(2)*z(3) + z(1)*x(2)*z(3);
option(prot);
option(redSB);option(redTail);option(mem);
ideal J = system("freegb",I,d,nvars(r));
lpDimCheck(J);
sickle(J,1,1,1,d); // dimension is 24872


LIB "fpadim.lib";
ring r = 0,(x,y,z),dp;
int d =10; // degree
def R = makeLetterplaceRing(d);
setring R;
ideal I = x(1)*y(2) + y(1)*z(2), x(1)*x(2) + x(1)*y(2) - y(1)*x(2) - y(1)*y(2);
option(prot);
option(redSB);option(redTail);option(mem);
ideal J = system("freegb",I,d,3);
lpDimCheck(J);
sickle(J,1,1,1,d);
 */


/*
   Example for computing GK dimension:
   returns a ring which contains an ideal I
   run gkDim(I) inside this ring and it should return 2n (the GK dimension
   of n-th Weyl algebra including evaluation operators).

   proc createWeylEx(int n, int d)
   "
   "
   {
   int baseringdef;
   if (defined(basering)) // if a basering is defined, it should be saved for later use
   {
   def save = basering;
   baseringdef = 1;
   }
   ring r = 0,(d(1..n),x(1..n),e(1..n)),dp;
   def R = makeLetterplaceRing(d);
   setring R;
   ideal I; int i,j;

   for (i = 1; i <= n; i++)
   {
   for (j = i+1; j<= n; j++)
   {
   I[size(I)+1] = lpMult(var(i),var(j));
   }
   }

   for (i = 1; i <= n; i++)
   {
   for (j = i+1; j<= n; j++)
   {
   I[size(I)+1] = lpMult(var(n+i),var(n+j));
   }
   }
   for (i = 1; i <= n; i++)
   {
   for (j = 1; j<= n; j++)
   {
   I[size(I)+1] = lpMult(var(i),var(n+j));
   }
   }
   for (i = 1; i <= n; i++)
   {
   for (j = 1; j<= n; j++)
   {
   I[size(I)+1] = lpMult(var(i),var(2*n+j));
   }
   }
   for (i = 1; i <= n; i++)
   {
   for (j = 1; j<= n; j++)
   {
   I[size(I)+1] = lpMult(var(2*n+i),var(n+j));
   }
   }
   for (i = 1; i <= n; i++)
   {
   for (j = 1; j<= n; j++)
   {
   I[size(I)+1] = lpMult(var(2*n+i),var(2*n+j));
   }
   }
   I = simplify(I,2+4);
   I = letplaceGBasis(I);
   export(I);
   if (baseringdef == 1) {setring save;}
   return(R);
   }

proc TestGKAuslander3()
{
  ring r = (0,q),(z,x,y),(dp(1),dp(2));
  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
  R; setring R; // sets basering to Letterplace ring
  ideal I;
  I = q*x(1)*y(2) - y(1)*x(2), z(1)*y(2) - y(1)*z(2), z(1)*x(2) - x(1)*z(2);
  I = letplaceGBasis(I);
  lpGkDim(I); // must be 3
  I = x(1)*y(2)*z(3) - y(1)*x(2), z(1)*y(2) - y(1)*z(2), z(1)*x(2) - x(1)*z(2);//gkDim = 2
  I = letplaceGBasis(I); // not finite BUT contains a poly in x,y only
  lpGkDim(I); // must be 4

  ring r = 0,(y,x,z),dp;
  def R = makeLetterplaceRing(10); // constructs a Letterplace ring
  R; setring R; // sets basering to Letterplace ring
  ideal I;
  I = x(1)*y(2)*z(3) - y(1)*x(2), z(1)*y(2) - y(1)*z(2), z(1)*x(2) - x(1)*z(2);//gkDim = 2
  I = letplaceGBasis(I); // computed as it would be homogenized; infinite
  poly p = x(1)*y(2)*y(3)*x(4)-y(1)*x(2)*x(3)*y(4);
  lpNF(p, I); // 0 as expected

  // with inverse of z
  ring r = 0,(iz,z,x,y),dp;
  def R = makeLetterplaceRing(11); // constructs a Letterplace ring
  R; setring R; // sets basering to Letterplace ring
  ideal I;
  I = x(1)*y(2)*z(3) - y(1)*x(2), z(1)*y(2) - y(1)*z(2), z(1)*x(2) - x(1)*z(2),
    iz(1)*y(2) - y(1)*iz(2), iz(1)*x(2) - x(1)*iz(2), iz(1)*z(2)-1, z(1)*iz(2) -1;
  I = letplaceGBasis(I); //
  setring r;
  def R2 = makeLetterplaceRing(23); // constructs a Letterplace ring
  setring R2; // sets basering to Letterplace ring
  ideal I = imap(R,I);
  lpGkDim(I);


  ring r = 0,(t,z,x,y),(dp(2),dp(2));
  def R = makeLetterplaceRing(20); // constructs a Letterplace ring
  R; setring R; // sets basering to Letterplace ring
  ideal I;
  I = x(1)*y(2)*z(3) - y(1)*x(2)*t(3), z(1)*y(2) - y(1)*z(2), z(1)*x(2) - x(1)*z(2),
    t(1)*y(2) - y(1)*t(2), t(1)*x(2) - x(1)*t(2), t(1)*z(2) - z(1)*t(2);//gkDim = 2
  I = letplaceGBasis(I); // computed as it would be homogenized; infinite
  LIB "elim.lib";
  ideal Inoz = nselect(I,intvec(2,6,10,14,18,22,26,30));
  for(int i=1; i<=20; i++)
  {
    Inoz=subst(Inoz,t(i),1);
  }
  ideal J = x(1)*y(2)*y(3)*x(4)-y(1)*x(2)*x(3)*y(4);
  J = letplaceGBasis(J);

  poly p = x(1)*y(2)*y(3)*x(4)-y(1)*x(2)*x(3)*y(4);
  lpNF(p, I); // 0 as expected

  ring r2 = 0,(x,y),dp;
  def R2 = makeLetterplaceRing(50); // constructs a Letterplace ring
  setring R2;
  ideal J = x(1)*y(2)*y(3)*x(4)-y(1)*x(2)*x(3)*y(4);
  J = letplaceGBasis(J);
}

*/


/*   actual work:
// downup algebra A
LIB "fpadim.lib";
ring r = (0,a,b,g),(x,y),Dp;
def R = makeLetterplaceRing(6); // constructs a Letterplace ring
setring R;
poly F1 = g*x(1);
poly F2 = g*y(1);
ideal J = x(1)*x(2)*y(3)-a*x(1)*y(2)*x(3) - b*y(1)*x(2)*x(3) - F1,
x(1)*y(2)*y(3)-a*y(1)*x(2)*y(3) - b*y(1)*y(2)*x(3) - F2;
J = letplaceGBasis(J);
lpGkDim(J); // 3 == correct

// downup algebra B
LIB "fpadim.lib";
ring r = (0,a,b,g, p(1..7),q(1..7)),(x,y),Dp;
def R = makeLetterplaceRing(6); // constructs a Letterplace ring
setring R;
ideal imn = 1, y(1)*y(2)*y(3), x(1)*y(2), y(1)*x(2), x(1)*x(2), y(1)*y(2), x(1), y(1);
int i;
poly F1, F2;
for(i=1;i<=7;i++)
{
F1 = F1 + p(i)*imn[i];
F2 = F2 + q(i)*imn[i];
}
ideal J = x(1)*x(2)*y(3)-a*x(1)*y(2)*x(3) - b*y(1)*x(2)*x(3) - F1,
x(1)*y(2)*y(3)-a*y(1)*x(2)*y(3) - b*y(1)*y(2)*x(3) - F2;
J = letplaceGBasis(J);
lpGkDim(J); // 3 == correct

 */