Changeset c8fd40 in git

Timestamp:

Nov 9, 2018, 4:09:43 PM (6 years ago)

Author:

Karim Abou Zeid <karim23697@…>

Branches:

(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd08f5f0bb3329b8ca19f23b74cb1473686415c3a')

Children:

6fac22d0e553731fefbeacc9e5e1cc9683b38ec2

Parents:

4aa50a694d58bf7ccd9553bb5644b3297c61b608

git-author:

Karim Abou Zeid <karim23697@gmail.com>2018-11-09 16:09:43+01:00

git-committer:

Karim Abou Zeid <karim23697@gmail.com>2018-11-14 16:38:23+01:00

Message:

Kernel stuff for uf graph

Files:

• Singular/LIB/fpaprops.lib (modified) (13 diffs)
• Singular/extra.cc (modified) (1 diff)
• libpolys/polys/shiftop.cc (modified) (1 diff)
• libpolys/polys/shiftop.h (modified) (1 diff)

Singular/LIB/fpaprops.lib

@@ -622 +622 @@
   // iterate through all vertices
 
-  intvec visited;
-  visited[n] = 0;
-
-  intvec cyclic;
-  cyclic[n] = 0;
+  intvec visited = 0:n;
+  intvec cyclic = 0:n;
+  intvec countedCycles = -2:n;
 
   int maxCycleCount = 0;
   for (int v = 1; v <= n; v++) {
-    int cycleCount = countCycles(UG, v, visited, cyclic, 0);
-    if (cycleCount == -1) {
+    countedCycles = countCycles(UG, v, visited, cyclic, 0, countedCycles);
+    dbprint("counted: " + string(v) + "/" + string(n));
+    if (countedCycles[v] == -1) {
       return(-1);
     }
-    if (cycleCount > maxCycleCount) {
-      maxCycleCount = cycleCount;
-    }
-    kill cycleCount;
+    if (countedCycles[v] > maxCycleCount) {
+      maxCycleCount = countedCycles[v];
+    }
   } kill v;
-  return(maxCycleCount);
-}
-
-static proc countCycles(intmat G, int v, intvec visited, intvec cyclic, intvec path)
+  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.
@@ -651 +649 @@
 "
 {
+  if (countedCycles[v] > -2) {
+    return (countedCycles);
+  }
   // Mark the current vertex as visited
   visited[v] = 1;
@@ -669 +670 @@
           if(cyclic[path[j]] == 1) {
             // 1.1 if yes, return -1
-            return (-1);
+            countedCycles[v] = -1;
+            return (countedCycles);
           }
           if (path[j] == w) {
@@ -694 +696 @@
         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);
+          countedCycles = countCycles(G, path[j], visited, cyclic, path, countedCycles);
+          if(countedCycles[path[j]] == -1) {
+            countedCycles[v] = -1;
+            return (countedCycles);
           }
-          if (cycleCount > maxCycleCount) {
-            maxCycleCount = cycleCount;
+          if (countedCycles[path[j]] > maxCycleCount) {
+            maxCycleCount = countedCycles[path[j]];
           }
-          kill cycleCount;
           if (path[j] == w) {
             break;
@@ -711 +713 @@
         kill maxCycleCount;
       } else {
-        int cycleCount = countCycles(G, w, visited, cyclic, path);
-        if (cycleCount == -1) {
-          return(-1);
-        }
-        if (cycleCount > cycles) {
-          cycles = cycleCount;
-        }
-        kill cycleCount;
+        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;
   // printf("Path: %s countCycles: %s", path, cycles); // DEBUG
-  return(cycles);
+  countedCycles[v] = cycles;
+  return (countedCycles);
 }
 
@@ -738 +741 @@
 "
 {
+  dbprint("computing maxDeg");
   int l = maxDeg(G);
-  list LG = lpId2ivLi(G);
-  list SW = ivStandardWords(LG, l - 1); // vertices
-  int n = size(SW);
+  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 = attrib(basering, "isLetterplaceRing");
+  // 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]
-      intvec v = SW[i];
-      intvec w = SW[j];
+      poly v = SW[i];
+      poly 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)) {
+        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 && !lpdivides(G, v * lpVarAt(w, l - 1))) {
         UG[i,j] = 1;
       }
-      kill v; kill w; kill v_overlap; kill w_overlap; kill vw;
+      kill v; kill w; kill v_overlap; kill w_overlap;
     } kill j;
   } kill i;
@@ -799 +807 @@
   } kill i;
   return (l);
+}
+
+static proc lpVarAt(poly p, int pos) {
+  return (system("lpVarAt", p, pos));
+}
+
+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 (lpdivides(G, words[i])) {
+      words = delete(words, i);
+    }
+  } kill i;
+  return (words);
 }
 
@@ -841 +870 @@
 }
 
+proc lpdivides(ideal I, poly p) {
+  for (int i = 1; i <= ncols(I); i++) {
+    if (system("lpdivides", I[i], p)) {
+      return (1);
+    }
+  } kill i;
+  return (0);
+}
+
 static proc ivdivides(list G, intvec iv) {
   for (int k = 1; k <= size(G); k++) {
@@ -858 +896 @@
 RETURN: int
 PURPOSE: Determines the Gelfand Kirillov dimension of A/<G>
-@*       -1 means it is positive infinite
+@*       -1 means positive infinite
 ASSUME: - basering is a Letterplace ring
 @*      - G is a Groebner basis
@@ -877 +915 @@
     // also if G is the whole ring return minus infinity
     if (leadmonom(G[i]) == 1) {
-      ERROR("Gk-Dim not defined for 0-ring")
+      ERROR("GK-Dim not defined for 0-ring")
     }
     kill d;
@@ -896 +934 @@
   }
 
+  dbprint("computing Ufnarovskii graph");
   intmat UG = lpUfGraph(G);
+  dbprint("Ufnarovskii graph: ", string(UG));
 
   // check special case 2
@@ -905 +945 @@
 
   // check special case 3
+  dbprint("topological sorting of Ufnarovskii graph");
   UG = topologicalSort(UG);
 
+  dbprint("check if Ufnarovskii graph is DAG");
   if (imIsUpRightTriangle(UG)) {
     UG = eliminateZerosUpTriangle(UG);
@@ -912 +954 @@
       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));
 }

Singular/extra.cc

@@ -1209 +1209 @@
     else
   #endif
+  /*==================== divide-test for freeGB  =================*/
+  #ifdef HAVE_SHIFTBBA
+    if (strcmp(sys_cmd, "lpdivides") == 0)
+    {
+      const short t[]={2,POLY_CMD,POLY_CMD};
+      if (iiCheckTypes(h,t,1))
+      {
+        poly p=(poly)h->CopyD();
+        poly q=(poly)h->next->CopyD();
+        res->rtyp = INT_CMD;
+        res->data = (void*)(long)p_LPDivisibleBy(p, q, currRing);
+        return FALSE;
+      }
+      else return TRUE;
+    }
+    else
+  #endif
+  /*==================== get var for freeGB  ====================*/
+  #ifdef HAVE_SHIFTBBA
+    if (strcmp(sys_cmd, "lpVarAt") == 0)
+    {
+      const short t[]={2,POLY_CMD,INT_CMD};
+      if (iiCheckTypes(h,t,1))
+      {
+        poly p=(poly)h->CopyD();
+        int pos=(int)((long)(h->next->Data()));
+        res->rtyp = POLY_CMD;
+        res->data = p_LPVarAt(p, pos, currRing);
+        return FALSE;
+      }
+      else return TRUE;
+    }
+    else
+  #endif
   /*==================== pcv ==================================*/
   #ifdef HAVE_PCV

libpolys/polys/shiftop.cc

@@ -659 +659 @@
   omFreeSize((ADDRESS) B, (b+1)*sizeof(int));
   return(1);
+}
+
+BOOLEAN p_LPDivisibleBy(poly a, poly b, const ring r)
+{
+  pIfThen1(b!=NULL, p_LmCheckPolyRing1(b, r));
+  pIfThen1(a!=NULL, p_LmCheckPolyRing1(a, r));
+
+  if (b == NULL) return TRUE;
+  if (a != NULL && (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r)))
+      return _p_LPLmDivisibleByNoComp(a,b,r);
+  return FALSE;
+}
+
+BOOLEAN p_LPLmDivisibleBy(poly a, poly b, const ring r)
+{
+  p_LmCheckPolyRing1(b, r);
+  pIfThen1(a != NULL, p_LmCheckPolyRing1(b, r));
+  if (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r))
+    return _p_LPLmDivisibleByNoComp(a, b, r);
+  return FALSE;
+}
+
+BOOLEAN _p_LPLmDivisibleByNoComp(poly a, poly b, const ring r)
+{
+  if(p_LmIsConstantComp(a, r))
+    return TRUE;
+#ifdef SHIFT_MULT_COMPAT_MODE
+  a = p_Head(a, r);
+  p_mLPunshift(a, r);
+  b = p_Head(b, r);
+  p_mLPunshift(b, r);
+#endif
+  int i = (r->N / r->isLPring) - p_LastVblock(a, r);
+  do {
+    int j = r->N - (i * r->isLPring);
+    bool divisible = true;
+    do
+    {
+      if (p_GetExp(a, j, r) > p_GetExp(b, j + (i * r->isLPring), r))
+      {
+        divisible = false;
+        break;
+      }
+      j--;
+    }
+    while (j);
+    if (divisible) return TRUE;
+    i--;
+  }
+  while (i > -1);
+#ifdef SHIFT_MULT_COMPAT_MODE
+  p_Delete(&a, r);
+  p_Delete(&b, r);
+#endif
+  return FALSE;
+}
+
+poly p_LPVarAt(poly p, int pos, const ring r)
+{
+  if (p == NULL || pos < 1 || pos > (r->N / r->isLPring)) return NULL;
+  poly v = p_One(r);
+  for (int i = (pos-1) * r->isLPring + 1; i <= pos * r->isLPring; i++) {
+    if (p_GetExp(p, i, r)) {
+      p_SetExp(v, i - (pos-1) * r->isLPring, 1, r);
+      return v;
+    }
+  }
+  return v;
 }
 

libpolys/polys/shiftop.h

@@ -48 +48 @@
 #define pmIsInV(p) p_mIsInV(p, currRing)
 
+BOOLEAN p_LPDivisibleBy(poly a, poly b, const ring r);
+BOOLEAN p_LPLmDivisibleBy(poly a, poly b, const ring r);
+BOOLEAN _p_LPLmDivisibleByNoComp(poly a, poly b, const ring r);
+
+poly p_LPVarAt(poly p, int pos, const ring r);
+
 ring freeAlgebra(ring r, int d);
 #endif