////////////////////////////////////////////////////////////////
version="version fpaprops.lib 4.1.2.0 Feb_2019 "; // $Id$
category="Noncommutative";
info="
LIBRARY: fpaprops.lib   Algorithmic ring-theoretic properties of finitely presented algebras (Letterplace)
AUTHORS: Karim Abou Zeid,       karim.abou.zeid at rwth-aachen.de

Support: Project II.6 in the transregional collaborative research centre
SFB-TRR 195 'Symbolic Tools in Mathematics and their Application' of the German DFG

OVERVIEW:
In this library, algorithms for computing various ring-theoretic properties of
finitely presented algebras are implemented.
Applicability: Letterplace rings.

REFERENCES:
Huishi Li: Groebner bases in ring theory. World Scientific, 2010.

KEYWORDS: finitely presented algebra; ring theory; Letterplace Groebner basis;
growth of algebra; Gelfand-Kirillov dimension; global homological dimension; semi-prime ideal; Ufnarovski graph

SEE ALSO:  fpadim_lib, freegb_lib

PROCEDURES:
  lpNoetherian(<GB>);     check whether A/<LM(GB)> is (left/right) Noetherian
  lpIsSemiPrime(<GB>);    check whether A/<LM(GB)> is semi-prime
  lpIsPrime(<GB>);        check whether A/<LM(GB)> is prime
  lpGkDim(<GB>);          alias for dim(<GB>)
  teach_lpGkDim(<GB>);    deprecated, kept for teaching purposes. use dim(<GB>) instead.
  lpGlDimBound(<GB>);     compute an upper bound for the global dimension of A/<GB>
  lpSubstitute();         substitute a variable with polynomials (ring homomorphism)
  lpCalcSubstDegBound();  utility for lpSubstitute
";

LIB "fpadim.lib";

/* very fast and cheap test of consistency and functionality
  DO NOT make it static !
  after adding the new proc, add it here */
proc tstfpaprops()
{
    example lpNoetherian;
    example lpIsSemiPrime;
    example lpIsPrime;
    example lpGlDimBound;
    example lpGkDim;
    example lpSubstitute;
    example lpCalcSubstDegBound;
};


////////////////////////////////////////////////////////////////////
proc lpNoetherian(ideal G)
"USAGE: lpNoetherian(G); G an ideal in a Letterplace ring
RETURN: int
      0 not Noetherian
      1 left Noetherian
      2 right Noetherian
      3 Noetherian
      4 weak Noetherian
PURPOSE: Check whether the monomial algebra A/<LM(G)> is (left/right) noetherian
ASSUME: - basering is a Letterplace ring
      - G is a Groebner basis
THEORY: lpNoetherian works with the monomial algebra A/<LM(G)>.
If it gives an affirmative answer for one of the properties, then it
holds for both A/<LM(G)> and A/<G>. However, a negative answer applies
only to A/<LM(G)> and not necessarily to A/<G>.
NOTE: Weak Noetherian means that two-sided ideals in A/<LM(G)> satisfy
the acc (ascending chain condition).
"
{
  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) {
      ERROR("Noetherianity not defined for 0-ring")
    }
    kill d;
  } kill i;
  // if longest word has length 1 we handle it as a special case
  if (l == 1) {
    int n = lpVarBlockSize(basering); // variable count
    int ncgenCount = lpNcgenCount(basering);
    int k = size(G);
    if (k == n - ncgenCount) { // only the field left
      return(3); // every field is noetherian
    }
    if (k == n - ncgenCount - 1) { // V = {1} with loop
      return(3);
    }
    if (k <= n ncgenCount - 2) { // V = {1} with more than one loop
      return(0);
    }
  }

  intmat UG = lpUfnarovskiGraph(G)[1];

  // 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 incoming/outcom
  intvec visited;
  visited[ncols(UG)] = 0;
  int inFlag, outFlag, inOutFlag;
  for (int v = 1; v <= ncols(UG) && (inFlag + outFlag) != 3; v++) {
    int inOutFlags = inOutCommingEdgesInCycles(UG, v, visited, 0);
    if (inOutFlags == 1) {
      inFlag = 1;
    }
    if (inOutFlags == 2) {
      outFlag = 1;
    }
    if (inOutFlags == 3) {
      inFlag = 1;
      outFlag = 1;
    }
    if (inOutFlags == 4) {
      inOutFlag = 1;
    }
    if (inOutFlags == 5) {
      inFlag = 1;
      inOutFlag = 1;
    }
    if (inOutFlags == 6) {
      outFlag = 1;
      inOutFlag = 1;
    }
    if (inOutFlags == 7) {
      inFlag = 1;
      outFlag = 1;
      inOutFlag = 1;
    }
    kill inOutFlags;
  } kill v;
  int noetherian = 3 - 1*inFlag - 2*outFlag;
  if (noetherian == 0) {
    return (4 - 4*inOutFlag); // weak noetherian
  }
  return (noetherian);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = freeAlgebra(r, 5);
  setring R;
  ideal G = x*x, y*x; // K<x,y>/<xx,yx> is right noetherian
  lpNoetherian(G);
}

static proc inOutCommingEdgesInCycles(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, inOutFlag;

  for (int w = 1; w <= ncols(G) && (inFlag + outFlag) != 3; w++) {
    if (G[v,w] == 1) {
      if (visited[w] == 1) { // new cycle
        int tmpInFlag;
        int tmpOutFlag;
        if (v == w) { // cycle is a loop
          for (int u = 1; u <= ncols(G); u++) {
            if (G[v,u] && u != v) {
              outFlag = 1;
              tmpOutFlag = 1;
            }
            if (G[u,v] && u != v) {
              inFlag = 1;
              tmpInFlag = 1;
            }
          } kill u;
        } 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 = 1;
                    tmpOutFlag = 1;
                  }
                } else {
                  if (u != w) {
                    outFlag = 1;
                    tmpOutFlag = 1;
                  }
                }
              }
              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;
                    tmpInFlag = 1;
                  }
                } else {
                  if (u != v) {
                    inFlag = 1;
                    tmpInFlag = 1;
                  }
                }
              }
            } kill u;
            if (path[i] == w) {
              break;
            }
          } kill i;
        }
        if (tmpInFlag > 0 && tmpOutFlag > 0) {
          // there are both in and outcomming edges in this cycle
          inOutFlag = 1;
        }
        kill tmpInFlag;
        kill tmpOutFlag;
      } else {
        int inOutFlags = inOutCommingEdgesInCycles(G, w, visited, path);
        if (inOutFlags == 1) {
          inFlag = 1;
        }
        if (inOutFlags == 2) {
          outFlag = 1;
        }
        if (inOutFlags == 3) {
          inFlag = 1;
          outFlag = 1;
        }
        if (inOutFlags == 4) {
          inOutFlag = 1;
        }
        if (inOutFlags == 5) {
          inFlag = 1;
          inOutFlag = 1;
        }
        if (inOutFlags == 6) {
          outFlag = 1;
          inOutFlag = 1;
        }
        if (inOutFlags == 7) {
          inFlag = 1;
          outFlag = 1;
          inOutFlag = 1;
        }
        kill inOutFlags;
      }
    }
  } kill w;

  return (1*inFlag + 2*outFlag + 4*inOutFlag);
}

proc lpIsSemiPrime(ideal G)
"USAGE: lpIsSemiPrime(G); G an ideal in a Letterplace ring
RETURN: boolean
PURPOSE: Check whether A/<LM(G)> is semi-prime ring,
alternatively whether <LM(G)> is a semi-prime ideal in A.
ASSUME: - basering is a Letterplace ring
      - G is a Groebner basis
THEORY: A (two-sided) ideal I in the ring A is semi-prime, if for any a in A one has
aAa subseteq I implies a in I.
NOTE: lpIsSemiPrime works with the monomial algebra A/<LM(G)>.
A positive answer holds for both A/<LM(G)> and A/<G>, while
a negative answer applies only to A/<LM(G)> and not necessarily to
A/<G>.
"
{
    // old theory part: that is when p * (A/<LM(G)>) * p != 0 for all p in (A/<LM(G)> - {0}).
  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) {
      ERROR("primeness not defined for 0-ring")
    }
    kill d;
  } kill i;
  // if longest word has length 1 we handle it as a special case
  if (l == 1) {
    return(1);
  }

  list VUG = lpUfnarovskiGraph(G);
  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);
    }
  } kill i;

  return (1);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x1,x2),dp;
  def R = freeAlgebra(r, 5);
  setring R;
  ideal G = x1*x2, x2*x1; // K<x1,x2>/<x1*x2,x2*x1> is semi prime
  lpIsSemiPrime(G);
}

// checks whether a monomial is a cyclic monomial
static 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);
        }
      }
    } kill i;
  } 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 A/<LM(G)> is prime ring,
alternatively whether <LM(G)> is a prime ideal in A.
ASSUME: - basering is a Letterplace ring
      - G is a Groebner basis
THEORY: A (two-sided) ideal I in the ring A is prime, if for any a,b in A one has
aAb subseteq I implies a in I or b in I.
NOTE: lpIsPrime works with the monomial algebra A/<LM(G)>.
A positive answer holds for both A/<LM(G)> and A/<G>, while
a negative answer applies only to A/<LM(G)> and not necessarily to A/<G>.
"
{
    // old theory part: that is when p1 * (A/<LM(G)>) * p2 != 0 for all p1, p2 in (A/<LM(G)> - {0}).
  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) {
      ERROR("primeness not defined for 0-ring")
    }
    kill d;
  } kill i;
  // if longest word has length 1 we handle it as a special case
  if (l == 1) {
    return(1);
  }

  list VUG = lpUfnarovskiGraph(G);
  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);
      }
    } kill j;
  } kill i;

  // 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;
      }
    } kill j;
  } kill i;

  return (1);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  def R = freeAlgebra(r, 5);
  setring R;
  ideal G = x*x, y*y; // 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);
      }
    }
  } kill i;

  return (0);
}

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
  } kill i;
  intmat UGk = UG;
  intmat zero[nrows(UGk)][ncols(UGk)];
  int k = 1;
  while (UGk != zero) {
    UGk = UGk * UG;
    k++;
  }
  return (k);
}

static 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); }
    } kill j;
  } kill i;
  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
            }
          } kill k;
        }
      } kill j;
      G = imDelRowCol(G,i,i); // remove vertex i
    }
  } kill 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];
        kill newi; kill newj;
      }
    } kill j;
  } kill i;
  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;
        }
      } kill k;
      if (incoming == 0) {
        v = j;
        kill incoming;
        break;
      } else {
        if (j == n) {
          // G contains at least one cycle, abort
          return (G);
        }
      }
      kill incoming;
    } kill j;

    // swap v and i
    if (v != i) {
      G = imPermcol(G, v, i);
      G = imPermrow(G, v, i);
    }
    kill v;
  } kill 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 = 0:n;
  intvec cyclic = 0:n;
  intvec countedCycles = -2:n;

  int maxCycleCount = 0;
  for (int v = 1; v <= n; v++) {
    countedCycles = countCycles(UG, v, visited, cyclic, 0, countedCycles);
    dbprint("counted " + string(countedCycles[v]) + " cycles from vertex " + string(v) + "/" + string(n) + " (cache: " + string(countedCycles) + ")");
    if (countedCycles[v] == -1) {
      return(-1);
    }
    if (countedCycles[v] > maxCycleCount) {
      maxCycleCount = countedCycles[v];
    }
  } kill v;
  return (maxCycleCount);
}

static proc countCycles(intmat G, int v, intvec visited, intvec cyclic, intvec path, intvec countedCycles)
"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
"
{
  if (countedCycles[v] > -2) {
    return (countedCycles);
  }
  // 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) { // found new cycle
        // 1. for all vertices in path until w, check if they are cyclic
        for (int j = size(path); j >= 1; j--) {
          if(cyclic[path[j]] == 1) {
            // 1.1 if yes, return -1
            countedCycles[v] = -1;
            return (countedCycles);
          }
          if (path[j] == w) {
            break;
          }
        } kill j;

        // 2. otherwise cycles++
        for (int j = size(path); j >= 1; j--) {
          // 2.2 remove the edges from that cycle and mark the vertices as cyclic
          if (j == size(path)) { // special case in the first 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;
          }
        } kill j;

        // 3. countCycles() on all these vertices
        int maxCycleCount = 0;
        for (int j = size(path); j >= 1; j--) {
          countedCycles = countCycles(G, path[j], visited, cyclic, path, countedCycles);
          if(countedCycles[path[j]] == -1) {
            countedCycles[v] = -1;
            return (countedCycles);
          }
          if (countedCycles[path[j]] > maxCycleCount) {
            maxCycleCount = countedCycles[path[j]];
          }
          if (path[j] == w) {
            break;
          }
        } kill j;
        if (maxCycleCount >= cycles) {
          cycles = maxCycleCount + 1;
        }
        kill maxCycleCount;
      } else {
        countedCycles = countCycles(G, w, visited, cyclic, path, countedCycles);
        if (countedCycles[w] == -1) {
          countedCycles[v] = -1;
          return (countedCycles);
        }
        if (countedCycles[w] > cycles) {
          cycles = countedCycles[w];
        }
      }
    }
  } kill w;
  countedCycles[v] = cycles;
  return (countedCycles);
}

// Ufnarovskii graph is now implemented in the dynamic module (freeAlgebra.cc)
/* proc lpUfnarovskiGraph(ideal G, list #) */
/* "USAGE: lpUfnarovskiGraph(G); G a set of monomials in a letterplace ring. */
/* RETURN: intmat or list */
/* NOTE: lpUfnarovskiGraph(G); returns intmat. lpUfnarovskiGraph(G,1); returns list L with L[1] an intmat and L[2] an ideal. */
/*       The intmat is the Ufnarovskij Graph and the ideal contains the vertices. */
/* 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 */
/* " */
/* { */
/*   dbprint("computing maxDeg"); */
/*   int l = maxDeg(G); */
/*   if (l - 1 == 0) { */
/*     // TODO: how should the graph look like when l - 1 = 0 ? */
/*     ERROR("Ufnarovskij graph not implemented for l = 1"); */
/*   } */
/*   int lV = lpVarBlockSize(basering); */
/*   // TODO: what if l <= 0? */
/*   dbprint("computing standard words"); */
/*   ideal SW = lpStandardWords(G, l - 1); // vertices */
/*   int n = ncols(SW); */
/*   dbprint("n = " + string(n)); */
/*   intmat UG[n][n]; // Ufnarovskij graph */
/*   for (int i = 1; i <= n; i++) { */
/*     for (int j = 1; j <= n; j++) { */
/*       dbprint("Ufnarovskii graph: " + string((i-1)*n + j) + "/" + string(n*n) + " entries"); */
/*       // [Studzinski page 76] */
/*       poly v = SW[i]; */
/*       poly w = SW[j]; */
/*       intvec v_overlap; */
/*       intvec w_overlap; */
/*       if (l - 1 > 1) { */
/*         v_overlap = leadexp(v); */
/*         w_overlap = leadexp(w); */
/*         v_overlap = v_overlap[(lV+1) .. (l-1)*lV]; */
/*         w_overlap = w_overlap[1 .. (l-2)*lV]; */
/*       } */
/*       if (v_overlap == w_overlap && !lpLmDivides(G, v * lpVarAt(w, l - 1))) { */
/*         UG[i,j] = 1; */
/*       } */
/*       kill v; kill w; kill v_overlap; kill w_overlap; */
/*     } kill j; */
/*   } kill i; */
/*   if (size(#) > 0) { */
/*     if (typeof(#[1]) == "int") { */
/*       if (#[1] != 0) { */
/*         list ret = UG; */
/*         ret[2] = SW; // the vertices */
/*         return (ret); */
/*       } */
/*     } */
/*   } */
/*   return (UG); */
/* } */
/* example */
/* { */
/*   "EXAMPLE:"; echo = 2; */
/*   ring r = 0,(x,y,z),dp; */
/*   def R = freeAlgebra(r, 5); // constructs a Letterplace ring */
/*   setring R; // sets basering to Letterplace ring */
/*   ideal I = x*y, x*z, z*y, z*z; */
/*   lpUfnarovskiGraph(I); */
/*   lpUfnarovskiGraph(I,1); */
/* } */

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;
    }
    kill d;
  } kill i;
  return (l);
}

static proc lpStandardWords(ideal G, int length)
"ASSUME: G is simplified
"
{
  if (length < 0) {
    return (delete(ideal(0), 1)); // no standard words
  }

  ideal words = maxideal(length);
  for (int i = ncols(words); i >= 1; i--) {
    if (lpLmDivides(G, words[i])) {
      words = delete(words, i);
    }
  } kill i;
  return (words);
}

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 nVars = lpVarBlockSize(basering) - lpNcgenCount(basering); // variable count
  list prevWords = ivStandardWords(G, length - 1);
  list words;
  for (int i = 1; i <= nVars; 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);
      }
      kill word;
    } kill j;
  } kill i;
  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);
  } kill 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);
      }
    }
  } kill k;
  return (0);
}

proc lpGkDim(ideal G)
"USAGE: lpGkDim(G); G an ideal in a letterplace ring
RETURN: int
PURPOSE: Determines the Gelfand Kirillov dimension of A/<G>
       -1 means positive infinite
ASSUME: - basering is a Letterplace ring
      - G is a Groebner basis
NOTE: Alias for dim(G)
"
{
  print("WARNING: `lpGkDim` is deprecated, you can use `dim` instead.");
  return (dim(G));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  ring R = freeAlgebra(r, 5);
  ideal I = z; // infinite GK dimension (-1)
  lpGkDim(I);
  I = x,y,z; I = std(I); // GK dimension 0
  lpGkDim(I);
  I = x*y, x*z, z*y, z*z; I = std(I); // GK dimension 2
  lpGkDim(I);
  ideal G = y*x - x*y, z*x - x*z, z*y - y*z; G = std(G);
  G;
  lpGkDim(G); // GK dimension 3
}

proc teach_lpGkDim(ideal G)
"USAGE: teach_lpGkDim(G); G an ideal in a letterplace ring
RETURN: int
PURPOSE: Determines the Gelfand Kirillov dimension of A/<G>
       -1 means positive infinite
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 return minus infinity
    if (leadmonom(G[i]) == 1) {
      ERROR("GK-Dim not defined for 0-ring")
    }
    kill d;
  } kill i;
  // if longest word has length 1, or G is the zero ideal, we handle it as a special case
  if (l == 1 || size(G) == 0) {
    int n = lpVarBlockSize(basering); // 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);
    }
  }

  dbprint("computing Ufnarovskii graph");
  intmat UG = lpUfnarovskiGraph(G)[1];
  if (printlevel >= voice + 1) {
    UG;
  }

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

  // check special case 3
  dbprint("topological sorting of Ufnarovskii graph");
  UG = topologicalSort(UG);
  if (printlevel >= voice + 1) {
    UG;
  }

  dbprint("check if Ufnarovskii graph is DAG");
  if (imIsUpRightTriangle(UG)) {
    UG = eliminateZerosUpTriangle(UG);
    if (ncols(UG) == 0 || nrows(UG) == 0) { // when the diagonal was zero
      return (0)
    }
    dbprint("DAG detected, using URTNZD growth alg");
    return(UfGraphURTNZDGrowth(UG));
  }

  // otherwise count cycles in the Ufnarovskij Graph
  dbprint("not a DAG, using regular growth alg");
  return(UfGraphGrowth(UG));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  def R = freeAlgebra(r, 5); // constructs a Letterplace ring
  R;
  setring R; // sets basering to Letterplace ring
  ideal I = z;//an example of infinite GK dimension
  teach_lpGkDim(I);
  I = x,y,z; // gkDim = 0
  teach_lpGkDim(I);
  I = x*y, x*z, z*y, z*z;//gkDim = 2
  teach_lpGkDim(I);
  ideal G = y*x - x*y, z*x - x*z, z*y - y*z; G = std(G);
  G;
  teach_lpGkDim(G); // 3, as expected for K[x,y,z]
}

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 Groebner 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]);
    }
  } kill 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);
  } kill i;
  G = simplify(G,2); // remove zero generators

  // delete variables in LM(G1) from the ring
  def save = basering;
  def R = basering;
  if (size(G1) > 0) {
    while (size(G1) > 0) {
      if (lpVarBlockSize(R) - lpNcgenCount(R) > 1) {
        def @R = R - string(G1[1]);
        R = @R;
        kill @R;
        setring R;
        /* ring R = lpDelVar(lp2iv(G1[1])[1]); // TODO replace with proper method */
        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
        G = twostd(G); // otherwise warning that G is no standard basis
        int gkDim = dim(G);
        if (gkDim >= 0) {
          return (gkDim);
        }
        kill gkDim;
      }
    }
  } kill i;

  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 = freeAlgebra(r, 5); // constructs a Letterplace ring
  setring R; // sets basering to Letterplace ring
  ideal G = x*x, y*y,x*y*x; // it is a monomial Groebner basis
  lpGlDimBound(G);
  ideal H = y*x - x*y; H = std(H); // H is a Groebner basis
  lpGlDimBound(H); // gl dim of K[x,y] is 2, as expected
}

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

static proc lpGraphOfNChains(ideal G) // G must be reduced
{
  list LG = lpId2ivLi(lead(G));
  int n = lpVarBlockSize(basering);
  int degbound = lpDegBound(basering);

  list V;
  for (int i = 0; i <= n; i++) {
    V[i+1] = i; // add 1 and all variables
  } kill i;
  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
      }
      kill v;
    } kill j;
    kill u;
  } kill i;
  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
  } kill i;

  // 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;
            }
          }
        } kill k;
      }
      kill uv;
    } kill j;
  } kill i;

  return(GNC);
}

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

proc lpSubstitute(poly f, ideal s1, ideal s2, list #)
"USAGE: lpSubstitute(f,s1,s2[,G]); f poly, 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
      - s1 contains a subset of the set of variables
      - s2 and s1 are of the same size
      - G is a Groebner basis,
      - the current ring has a sufficient degbound (which also can be calculated with lpCalcSubstDegBound())
NOTE: the procedure implements the image of a polynomial f
under an endomorphism of a free algebra, defined by s1 and s2:
variables, not present in s1, are left unchanged;
variable s1[k] is mapped to a polynomial s2[k].
- An optional ideal G extends the endomorphism as above to the morphism into the factor algebra K<X>/G.
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];
      if (varindex > 0) {
        int subindex = lpIndexOf(s1, var(varindex));
        if (subindex > 0) {
          s2[subindex] = lpNF(s2[subindex],G);
          fis = fis * s2[subindex];
        } else {
          fis = fis * lpNF(iv2lp(varindex),G);
        }
        /*fis = lpNF(fis,G);*/
        kill subindex;
      }
      kill varindex;
    } kill j;
    kill ivfi;
    fs = fs + fis;
    kill fis;
  }
  kill i;
  fs = lpNF(fs, G);
  return (fs);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  def R = freeAlgebra(r, 4);
  setring R;
  ideal G = x*y; // optional
  poly f = 3*x*x+y*x;
  ideal s1 = x, y;
  ideal s2 = y*z*z, x; // i.e. x --> yzz and y --> x
  // the substitution probably needs a higher degbound
  int minDegBound = lpCalcSubstDegBound(f,s1,s2);
  minDegBound; // thus the bound needs to be increased
  setring r; // back to original r
  def R1 = freeAlgebra(r, minDegBound);
  setring R1;
  lpSubstitute(imap(R,f), imap(R,s1), imap(R,s2));
  // the last parameter is optional; above it was G=<xy>
  // the output will be reduced with respect to G
  lpSubstitute(imap(R,f), imap(R,s1), imap(R,s2), imap(R,G));
}

// another example:
/*
  //////// EXAMPLE B ////////
  ring r = 0,(x,y,z),dp;
  def R = freeAlgebra(r, 4);
  setring R;

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

  int minDegBound = lpCalcSubstDegBound(I,s1,s2);
  setring r;
  def R1 = freeAlgebra(r, 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);
    }
  } kill i;
  return (-1);
}

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


static proc lpCalcSubstDegBoundSingle(poly f, ideal s1, ideal s2)
"USAGE: lpCalcSubstDegBoundSingle(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 lpCalcSubstDegBoundSingle; 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];
      if (varindex > 0) {
        int subindex = lpIndexOf(s1, var(varindex));
        if (subindex > 0) {
          tmpDegBound = tmpDegBound + deg(s2[subindex]);
        } else {
          tmpDegBound = tmpDegBound + 1;
        }
        kill subindex;
      }
      kill varindex;
    } kill j;
    if (tmpDegBound > maxDegBound) {
      maxDegBound = tmpDegBound;
    }
    kill ivfi; kill tmpDegBound;
  } kill i;

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

  return (maxDegBound);
}
example
{
  // see lpCalcSubstDegBound()
}

proc lpCalcSubstDegBound(ideal I, ideal s1, ideal s2)
"USAGE: lpCalcSubstDegBound(I,s1,s2); I ideal of polynomials, s1 ideal of variables to replace, s2 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 lpCalcSubstDegBound; shows examples
NOTE: convenience method
"
{
  int maxDegBound = 0;
  for (int i = 1; i <= size(I); i++) {
    int tmpDegBound = lpCalcSubstDegBoundSingle(I[i], s1, s2, #);
    if (tmpDegBound > maxDegBound) {
      maxDegBound = tmpDegBound;
    }
    kill tmpDegBound;
  } kill i;
  return (maxDegBound);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  def R = freeAlgebra(r, 4);
  setring R;
  ideal I = 3*x*x+y*x, x*y*x - z;
  ideal s1 = x, y; // z --> z
  ideal s2 = y*z*z, x; // i.e. x --> yzz and y --> x
  // the substitution probably needs a higher degbound
  lpCalcSubstDegBound(I,s1,s2);
  lpCalcSubstDegBound(I[1],s1,s2);
}

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);
    }
    kill IIF;
  } kill i;
  return(0);
}

// no longer working with new interface and new orderings
/* // TODO: use original ring attrib to create a new letterplace ring */
/* // removes a variable from a letterplace ring (a bit of a hack) */
/* static proc lpDelVar(int index) { */
/*   int lV = lpVarBlockSize(basering); // number of variables in the main block */
/*   int d = lpDegBound(basering); // degree bround */
/*   list LR = ringlist(basering); */

/*   if (!(index >= 1 && index <= lV)) { return (basering); } // invalid index */

/*   // remove from the variable list */
/*   for (int i = (d-1)*lV + index; i >= 1; i = i - lV) { */
/*     LR[2] = delete(LR[2], i); */
/*   } kill i; */

/*   // remove from a ordering */
/*   intvec aiv = LR[3][1][2]; */
/*   aiv = aiv[1..(d*lV-d)]; */
/*   LR[3][1][2] = aiv; */

/*   // remove block orderings */
/*   int del = (lV - index); */
/*   int cnt = -1; */
/*   for (int i = size(LR[3]); i >= 2; i--) { */
/*     if (LR[3][i][2] != 0) { */
/*       for (int j = size(LR[3][i][2]); j >= 1; j--) { */
/*         cnt++; // next 1 */
/*         if (cnt%lV == del) { */
/*           // delete */
/*           if (size(LR[3][i][2]) > 1) { // if we have more than one element left, delete one */
/*             LR[3][i][2] = delete(LR[3][i][2],j); */
/*           } else { // otherwise delete the whole block */
/*             LR[3] = delete(LR[3], i); */
/*             break; */
/*           } */
/*         } */
/*       } kill j; */
/*     } */
/*   } kill i; */

/*   def R = setLetterplaceAttributes(ring(LR),d,lV-1); */
/*   return (R); */
/* } */
/* example */
/* { */
/*   "EXAMPLE:"; echo = 2; */
/*   ring r = 0,(x,y,z),dp; */
/*   def A = freeAlgebra(r, 3); */
/*   setring A; A; */
/*   def R = lpDelVar(2); setring R; R; */
/* } */