source: git/Singular/LIB/fpaprops.lib @ 4d9f75e

spielwiese
Last change on this file since 4d9f75e was 4d9f75e, checked in by Karim Abou Zeid <karim23697@…>, 6 years ago
Fix console spam
  • Property mode set to 100644
File size: 31.7 KB
Line 
1////////////////////////////////////////////////////////////////
2version="version fpaprops.lib 1.0.0.0 Dec_2017 "; // $Id$
3category="Noncommutative";
4info="
5LIBRARY: fpaprops.lib   Algorithms for the properties of quotient algebras in the letterplace case
6AUTHORS: Karim Abou Zeid,       karim.abou.zeid at rwth-aachen.de
7
8Support: Project II.6 in the transregional collaborative research centre
9SFB-TRR 195 'Symbolic Tools in Mathematics and their Application' of the German DFG
10
11OVERVIEW:
12Algorithms for computing various properties of quotient algebras in the letterplace case.
13
14REFERENCES:
15Huishi Li: Groebner bases in ring theory. World Scientific, 2010.
16
17SEE ALSO:  fpadim_lib, freegb_lib
18
19PROCEDURES:
20  lpNoetherian(<GB>);     check whether A/<GB> is (left/right) noetherian
21  lpIsSemiPrime(<GB>);    check whether A/<GB> is semi prime
22  lpIsPrime(<GB>);        check whether A/<GB> is prime
23  lpGkDim(<GB>);          compute the Gelfand Kirillov dimension of A/<GB>
24  lpGlDimBound(<GB>);     compute an upper bound for the global dimension of A/<GB>
25  lpSubstitute();         substitute variable with polynomials
26  lpCalcSubstDegBound();  utility for lpSubstitute
27  lpCalcSubstDegBounds(); utility for lpSubstitute
28";
29
30LIB "fpadim.lib";
31////////////////////////////////////////////////////////////////////
32proc lpNoetherian(ideal G)
33"USAGE: lpNoetherian(G); G an ideal in a Letterplace ring
34RETURN: int
35@*      0 not noetherian
36@*      1 left noetherian
37@*      2 right noetherian
38@*      3 noetherian
39PURPOSE: Check whether R/<G> is (left/right) noetherian, where R is the basering
40ASSUME: - basering is a Letterplace ring
41@*      - G is a Groebner basis
42"
43{
44  G = lead(G);
45  G = simplify(G, 2+4+8);
46
47  // check special case 1
48  int l = 0;
49  for (int i = 1; i <= size(G); i++) {
50    // find the max degree in G
51    int d = deg(G[i]);
52    if (d > l) {
53      l = d;
54    }
55
56    // also if G is the whole ring
57    if (leadmonom(G[i]) == 1) {
58      ERROR("noetherianity not defined for 0-ring")
59    }
60    kill d;
61  } kill i;
62  // if longest word has length 1 we handle it as a special case
63  if (l == 1) {
64    int n = attrib(basering, "lV"); // variable count
65    int k = size(G);
66    if (k == n) { // only the field left
67      return(3); // every field is noetherian
68    }
69    if (k == n-1) { // V = {1} with loop
70      return(3);
71    }
72    if (k <= n-2) { // V = {1} with more than one loop
73      return(0);
74    }
75  }
76
77  intmat UG = lpUfGraph(G);
78
79  // check special case 2
80  intmat zero[nrows(UG)][ncols(UG)];
81  if (UG == zero) {
82    return (3);
83  }
84
85  if (!imHasLoops(UG) && imIsUpRightTriangle(topologicalSort(UG))) {
86    // UG is a DAG
87    return (3);
88  }
89
90  // DFS from every vertex, if cycle is found, check every vertex for incomming/outcom
91  intvec visited;
92  visited[ncols(UG)] = 0;
93  int inFlag, outFlag;
94  for (int v = 1; v <= ncols(UG) && (inFlag + outFlag) != 3; v++) {
95    int inOutFlags = inOrOutCommingEdgeInCycle(UG, v, visited, 0);
96    if (inOutFlags == 1) {
97      inFlag = 1;
98    }
99    if (inOutFlags == 2) {
100      outFlag = 2;
101    }
102    if (inOutFlags == 3) {
103      inFlag = 1;
104      outFlag = 2;
105    }
106    kill inOutFlags;
107  } kill v;
108  return (3 - inFlag - outFlag);
109}
110example {
111  "EXAMPLE:"; echo = 2;
112  ring r = 0,(x,y),dp;
113  def R = makeLetterplaceRing(5);
114  setring R;
115  ideal G = x(1)*x(2), y(1)*x(2); // K<x,y>/<xx,yx> is right noetherian
116  lpNoetherian(G);
117}
118
119static proc inOrOutCommingEdgeInCycle(intmat G, int v, intvec visited, intvec path) {
120  // Mark the current vertex as visited
121  visited[v] = 1;
122
123  // Store the current vertex in path
124  if (path[1] == 0) {
125    path[1] = v;
126  } else {
127    path[size(path) + 1] = v;
128  }
129
130  int inFlag, outFlag;
131
132  for (int w = 1; w <= ncols(G) && (inFlag + outFlag) != 3; w++) {
133    if (G[v,w] == 1) {
134      if (visited[w] == 1) {
135        // new cycle
136        if (v == w) {
137          for (int u = 1; u <= ncols(G); u++) {
138            if (G[v,u] && u != v) {
139              outFlag = 2;
140            }
141            if (G[u,v] && u != v) {
142              inFlag = 1;
143            }
144          } kill u;
145        } else {
146          for (int i = size(path); i >= 1; i--) { // for each vertex in the path
147            // check for neighbors not directly next or prev in cycle
148            for (int u = 1; u <= ncols(G); u++) {
149              if (G[path[i],u] == 1) { // there is an edge to u
150                if (path[i] != v) {
151                  if (u != path[i+1]) { // and u is not the next element in the cycle
152                    outFlag = 2;
153                  }
154                } else {
155                  if (u != w) {
156                    outFlag = 2;
157                  }
158                }
159              }
160              if (G[u,path[i]] == 1) { // there is an edge from u
161                if (path[i] != w) {
162                  if (u != path[i-1]) { // and u is not the previous element in the cylce
163                    inFlag = 1;
164                  }
165                } else {
166                  if (u != v) {
167                    inFlag = 1;
168                  }
169                }
170              }
171            } kill u;
172            if (path[i] == w) {
173              break;
174            }
175          } kill i;
176        }
177      } else {
178        int inOutFlags = inOrOutCommingEdgeInCycle(G, w, visited, path);
179        if (inOutFlags == 1) {
180          inFlag = 1;
181        }
182        if (inOutFlags == 2) {
183          outFlag = 2;
184        }
185        if (inOutFlags == 3) {
186          inFlag = 1;
187          outFlag = 2;
188        }
189        kill inOutFlags;
190      }
191    }
192  } kill w;
193
194  return (inFlag + outFlag);
195}
196
197proc lpIsSemiPrime(ideal G)
198"USAGE: lpIsSemiPrime(G); G an ideal in a Letterplace ring
199RETURN: boolean
200PURPOSE: Check whether R/<G> is semi prime, where R is the basering
201ASSUME: - basering is a Letterplace ring
202      - G is a Groebner basis
203"
204{
205  G = lead(G);
206  G = simplify(G, 2+4+8);
207
208  // check special case 1
209  int l = 0;
210  for (int i = 1; i <= size(G); i++) {
211    // find the max degree in G
212    int d = deg(G[i]);
213    if (d > l) {
214      l = d;
215    }
216
217    // also if G is the whole ring
218    if (leadmonom(G[i]) == 1) {
219      ERROR("primeness not defined for 0-ring")
220    }
221    kill d;
222  } kill i;
223  // if longest word has length 1 we handle it as a special case
224  if (l == 1) {
225    return(1);
226  }
227
228  list VUG = lpUfGraph(G, 1);
229  intmat UG = VUG[1]; // the Ufnarovskij graph
230  ideal V = VUG[2]; // the vertices of UG (standard words with length = l-1)
231
232  list LG = lpId2ivLi(G);
233  list SW = ivStandardWordsUpToLength(LG, maxDeg(G));
234  list LV = lpId2ivLi(V);
235
236  // delete the 0 in SW
237  int indexofzero = ivIndexOf(SW, 0);
238  if (indexofzero > 0) { // should be always true when |SW| > 0
239    SW = delete(SW, indexofzero);
240  }
241
242  // check if each monomial in SW is cyclic
243  for (int i = 1; i <= size(SW); i++) {
244    if (!isCyclicInUfGraph(UG, LV, SW[i])) {
245      return (0);
246    }
247  } kill i;
248
249  return (1);
250}
251example {
252  "EXAMPLE:"; echo = 2;
253  ring r = 0,(x1,x2),dp;
254  def R = makeLetterplaceRing(5);
255  setring R;
256  ideal G = x1(1)*x2(2), x2(1)*x1(2); // K<x1,x2>/<x1*x2,x2*x1> is semi prime
257  lpIsSemiPrime(G);
258}
259
260// checks whether a monomial is a cyclic monomial
261static proc isCyclicInUfGraph(intmat UG, list LV, intvec u)
262{
263  if (ncols(UG) == 0) {return (0);} // UG is empty
264  if (u == 0) {return (0);} // 0 is never cyclic
265
266  int l = size(LV[1]) + 1;
267
268  int s = size(u);
269  if (s <= l - 1) {
270    for (int i = 1; i <= size(LV); i++) {
271      // for all vertices where u is a suffix
272      if(isSF(u, LV[i])) {
273        if (existsRoute(UG, i, i)) {
274          return (1);
275        }
276      }
277    } kill i;
278  } else { // size(u) > l - 1
279    int m = s - l + 1;
280
281    // there must be a route from v0 to vm
282    intvec v0 = u[1..(l-1)]; // first in route of u
283    intvec vm = u[m+1..m+(l-1)]; // last in route of u
284
285    int iv0 = ivIndexOf(LV, v0);
286    int ivm = ivIndexOf(LV, vm);
287    if (iv0 <= 0 || ivm <= 0) {
288      ERROR("u is not a standard word");
289    }
290
291    return (existsRoute(UG, ivm, iv0));
292  }
293
294  return (0);
295}
296
297proc lpIsPrime(ideal G)
298"USAGE: lpIsPrime(G); G an ideal in a Letterplace ring
299RETURN: boolean
300PURPOSE: Check whether R/<G> is prime, where R is the basering
301ASSUME: - basering is a Letterplace ring
302      - G is a Groebner basis
303"
304{
305  G = lead(G);
306  G = simplify(G, 2+4+8);
307
308  // check special case 1
309  int l = 0;
310  for (int i = 1; i <= size(G); i++) {
311    // find the max degree in G
312    int d = deg(G[i]);
313    if (d > l) {
314      l = d;
315    }
316
317    // also if G is the whole ring
318    if (leadmonom(G[i]) == 1) {
319      ERROR("primeness not defined for 0-ring")
320    }
321    kill d;
322  } kill i;
323  // if longest word has length 1 we handle it as a special case
324  if (l == 1) {
325    return(1);
326  }
327
328  list VUG = lpUfGraph(G, 1);
329  intmat UG = VUG[1]; // the Ufnarovskij graph
330  ideal V = VUG[2]; // the vertices of UG (standard words with length = l-1)
331
332  list LG = lpId2ivLi(G);
333  list LV = lpId2ivLi(V);
334
335  int n = ncols(UG);
336
337  // 1) for each vi vj there exists a route from vi to vj (means UG is connected)
338  for (int i = 1; i <= n; i++) {
339    for (int j = 1; j <= n; j++) {
340      if (!existsRoute(UG, i, j)) {
341        return (0);
342      }
343    } kill j;
344  } kill i;
345
346  // 2) any standard word with length < l-1 is a suffix of a vertex
347  list SW = ivStandardWordsUpToLength(LG, maxDeg(G) - 2); // < maxDeg - 1
348  if (size(SW) > 0 && size(LV) == 0) {return (0);}
349  for (int i = 1; i <= size(SW); i++) {
350    // check if SW[i] is a suffix of some LV
351    for (int j = 1; j <= size(LV); j++) {
352      if (!isSF(SW[i], LV[j])) {
353        if (j == size(LV)) {
354          return (0);
355        }
356      } else {
357        break;
358      }
359    }
360  } kill i;
361
362  return (1);
363}
364example {
365  "EXAMPLE:"; echo = 2;
366  ring r = 0,(x,y),dp;
367  def R = makeLetterplaceRing(5);
368  setring R;
369  ideal G = x(1)*x(2), y(1)*y(2); // K<x,y>/<xx,yy> is prime
370  lpIsPrime(G);
371}
372
373static proc existsRoute(intmat G, int v, int u, list #)
374"USAGE: existsRoute(G,v,u); G a graph, v and u vertices
375NOTE: don't pass anything to # (internal use for recursion)
376routes always have at least one edge
377"
378{
379  int n = ncols(G);
380
381  // init visited
382  intvec visited;
383  if (size(#) > 0) {
384    if (v == u) {return (1);} // don't check on first call so |route| >= 1 holds
385    visited = #[1];
386  } else { // first call
387    visited[n] = 0;
388  }
389
390  // mark current vertex as visited
391  visited[v] = 1;
392
393  // recursive DFS
394  for (int i = 1; i <= n; i++) {
395    if (G[v,i] && (!visited[i] || i == u)) { // i == u to allow routes from u to u
396      if (existsRoute(G, i, u, visited)) {
397        return (1);
398      }
399    }
400  } kill i;
401
402  return (0);
403}
404
405static proc UfGraphURTNZDGrowth(intmat UG) {
406  // URTNZD = upper right triangle non zero diagonal
407  for (int i = 1; i <= ncols(UG); i++) {
408    UG[i,i] = 0; // remove all loops
409  } kill i;
410  intmat UGk = UG;
411  intmat zero[nrows(UGk)][ncols(UGk)];
412  int k = 1;
413  while (UGk != zero) {
414    UGk = UGk * UG;
415    k++;
416  }
417  return (k);
418}
419
420static proc imIsUpRightTriangle(intmat M) {
421  for (int i = 1; i <= nrows(M); i++) {
422    for (int j = 1; j < i; j++) {
423      if(M[i,j] != 0) { return (0); }
424    } kill j;
425  } kill i;
426  return (1);
427}
428
429static proc eliminateZerosUpTriangle(intmat G) {
430  // G is expected to be an upper triangle matrix
431  for (int i = ncols(G); i >= 1; i--) { // loop order is important because we delete entries
432    if (G[i,i] == 0) { // i doesn't have a cycle
433      for (int j = 1; j < i; j++) {
434        if (G[j,i] == 1) { // j has an edge to i
435          for (int k = i + 1; k <= nrows(G); k++) {
436            if (G[i,k] == 1) {
437              G[j,k] = G[i,k]; // give j all edges from i
438            }
439          } kill k;
440        }
441      } kill j;
442      G = imDelRowCol(G,i,i); // remove vertex i
443    }
444  } kill i;
445  return (G);
446}
447
448static proc imDelRowCol(intmat M, int row, int col) {
449  // row and col are expected to be > 0
450  int nr = nrows(M);
451  int nc = ncols(M);
452  intmat Mdel[nr - 1][nc - 1];
453  for (int i = 1; i <= nr; i++) {
454    for (int j = 1; j <= nc; j++) {
455      if(i != row && j != col) {
456        int newi = i;
457        int newj = j;
458        if (i > row) { newi = i - 1; }
459        if (j > col) { newj = j - 1; }
460        Mdel[newi,newj] = M[i,j];
461        kill newi; kill newj;
462      }
463    } kill j;
464  } kill i;
465  return (Mdel);
466}
467
468static proc topologicalSort(intmat G) {
469  // NOTE: ignores loops
470  // NOTE: this takes O(|V^3|), can be optimized
471  int n = ncols(G);
472  for (int i = 1; i <= n; i++) { // only use the submat at i
473    // find a vertex v in the submat at i with no incoming edges
474    int v;
475    for (int j = i; j <= n; j++) {
476      int incoming = 0;
477      for (int k = i; k <= n; k++) {
478        if (k != j && G[k,j] == 1) {
479          incoming = 1;
480        }
481      } kill k;
482      if (incoming == 0) {
483        v = j;
484        break;
485      } else {
486        if (j == n) {
487          // G contains at least one cycle, abort
488          return (G);
489        }
490      }
491      kill incoming;
492    } kill j;
493
494    // swap v and i
495    if (v != i) {
496      G = imPermcol(G, v, i);
497      G = imPermrow(G, v, i);
498    }
499    kill v;
500  } kill i;
501  return (G);
502}
503
504static proc imPermcol (intmat A, int c1, int c2)
505{
506  intmat B = A;
507  int k = nrows(B);
508  B[1..k,c1] = A[1..k,c2];
509  B[1..k,c2] = A[1..k,c1];
510  return (B);
511}
512
513static proc imPermrow (intmat A, int r1, int r2)
514{
515  intmat B = A;
516  int k = ncols(B);
517  B[r1,1..k] = A[r2,1..k];
518  B[r2,1..k] = A[r1,1..k];
519  return (B);
520}
521
522static proc UfGraphGrowth(intmat UG)
523{
524  int n = ncols(UG); // number of vertices
525  // iterate through all vertices
526
527  intvec visited;
528  visited[n] = 0;
529
530  intvec cyclic;
531  cyclic[n] = 0;
532
533  int maxCycleCount = 0;
534  for (int v = 1; v <= n; v++) {
535    int cycleCount = countCycles(UG, v, visited, cyclic, 0);
536    if (cycleCount == -1) {
537      return(-1);
538    }
539    if (cycleCount > maxCycleCount) {
540      maxCycleCount = cycleCount;
541    }
542    kill cycleCount;
543  } kill v;
544  return(maxCycleCount);
545}
546
547static proc countCycles(intmat G, int v, intvec visited, intvec cyclic, intvec path)
548"USAGE: countCycles(G, v, visited, cyclic, path); G a Graph, v the vertex to
549start. The parameter visited, cyclic and path should be 0.
550RETURN: int
551@*:     Maximal number of distinct cycles
552PURPOSE: Calculate the maximal number of distinct cycles in a single path starting at v
553ASSUME: Basering is a Letterplace ring
554"
555{
556  // Mark the current vertex as visited
557  visited[v] = 1;
558
559  // Store the current vertex in path
560  if (path[1] == 0) {
561    path[1] = v;
562  } else {
563    path[size(path) + 1] = v;
564  }
565
566  int cycles = 0;
567  for (int w = 1; w <= ncols(G); w++) {
568    if (G[v,w] == 1) {
569      if (visited[w] == 1) { // found new cycle
570        // 1. for all vertices in path until w, check if they are cyclic
571        for (int j = size(path); j >= 1; j--) {
572          if(cyclic[path[j]] == 1) {
573            // 1.1 if yes, return -1
574            return (-1);
575          }
576          if (path[j] == w) {
577            break;
578          }
579        } kill j;
580
581        // 2. otherwise cycles++
582        for (int j = size(path); j >= 1; j--) {
583          // 2.2 remove the edges from that cycle and mark the vertices as cyclic
584          if (j == size(path)) { // special case in the first iteration
585            cyclic[v] = 1;
586            G[v, w] = 0;
587          } else {
588            cyclic[path[j]] = 1;
589            G[path[j], path[j+1]] = 0;
590          }
591          if (path[j] == w) {
592            break;
593          }
594        } kill j;
595
596        // 3. countCycles() on all these vertices
597        int maxCycleCount = 0;
598        for (int j = size(path); j >= 1; j--) {
599          int cycleCount = countCycles(G, path[j], visited, cyclic, path);
600          if(cycleCount == -1) {
601            return (-1);
602          }
603          if (cycleCount > maxCycleCount) {
604            maxCycleCount = cycleCount;
605          }
606          if (path[j] == w) {
607            break;
608          }
609          kill cycleCount;
610        } kill j;
611        if (maxCycleCount >= cycles) {
612          cycles = maxCycleCount + 1;
613        }
614        kill maxCycleCount;
615      } else {
616        int cycleCount = countCycles(G, w, visited, cyclic, path);
617        if (cycleCount == -1) {
618          return(-1);
619        }
620        if (cycleCount > cycles) {
621          cycles = cycleCount;
622        }
623        kill cycleCount;
624      }
625    }
626  } kill w;
627  // printf("Path: %s countCycles: %s", path, cycles); // DEBUG
628  return(cycles);
629}
630
631static proc lpUfGraph(ideal G, list #)
632"USAGE: lpUfGraph(G); G a set of monomials in a letterplace ring, # an optional parameter to return the vertex list when set
633RETURN: intmat
634PURPOSE: Constructs the Ufnarovskij graph induced by G
635      the adjacency matrix of the Ufnarovskij graph induced by G
636ASSUME: - basering is a Letterplace ring
637      - G are the leading monomials of a Groebner basis
638"
639{
640  int l = maxDeg(G);
641  list LG = lpId2ivLi(G);
642  list SW = ivStandardWords(LG, l - 1); // vertices
643  int n = size(SW);
644  intmat UG[n][n]; // Ufnarovskij graph
645  for (int i = 1; i <= n; i++) {
646    for (int j = 1; j <= n; j++) {
647      // [Studzinski page 76]
648      intvec v = SW[i];
649      intvec w = SW[j];
650      intvec v_overlap;
651      intvec w_overlap;
652      //TODO how should the graph look like when l - 1 = 0 ?
653      if (l - 1 == 0) {
654        ERROR("Ufnarovskij graph not implemented for l = 1");
655      }
656      if (l - 1 > 1) {
657        v_overlap = v[2 .. l-1];
658        w_overlap = w[1 .. l-2];
659      }
660      intvec vw = v;
661      vw[l] = w[l-1];
662      if (v_overlap == w_overlap && !ivdivides(LG, vw)) {
663        UG[i,j] = 1;
664      }
665      kill v; kill w; kill v_overlap; kill w_overlap; kill vw;
666    } kill j;
667  } kill i;
668  if (size(#) > 0) {
669    if (typeof(#[1]) == "int") {
670      if (#[1] == 1) {
671        list ret = UG;
672        ret[2] = ivL2lpI(SW); // the vertices
673        return (ret);
674      }
675    }
676  }
677  return (UG);
678}
679
680static proc maxDeg(ideal G)
681{
682  int l = 0;
683  for (int i = 1; i <= size(G); i++) { // find the max degree in G
684    int d = deg(G[i]);
685    if (d > l) {
686      l = d;
687    }
688    kill d;
689  } kill i;
690  return (l);
691}
692
693static proc ivStandardWords(list G, int length)
694"ASSUME: G is simplified
695"
696{
697  if (length <= 0) {
698    list words;
699    if (length == 0 && !ivdivides(G,0)) {
700      words[1] = 0; // iv = 0 means monom = 1
701    }
702    return (words); // no standard words
703  }
704  int lV = attrib(basering, "lV"); // variable count
705  list prevWords = ivStandardWords(G, length - 1);
706  list words;
707  for (int i = 1; i <= lV; i++) {
708    for (int j = 1; j <= size(prevWords); j++) {
709      intvec word = prevWords[j];
710      word[length] = i;
711      // assumes that G is simplified!
712      if (!ivdivides(G, word)) {
713        words = insert(words, word);
714      }
715      kill word;
716    } kill j;
717  } kill i;
718  return (words);
719}
720
721static proc ivStandardWordsUpToLength(list G, int length)
722"ASSUME: G is simplified
723"
724{
725  list words = ivStandardWords(G,0);
726  if (size(words) == 0) {return (words)}
727  for (int i = 1; i <= length; i++) {
728    words = words + ivStandardWords(G, i);
729  } kill i;
730  return (words);
731}
732
733static proc ivdivides(list G, intvec iv) {
734  for (int k = 1; k <= size(G); k++) {
735    if (isIF(G[k], iv)) {
736      return (1);
737    } else {
738      if (k == size(G)) {
739        return (0);
740      }
741    }
742  } kill k;
743  return (0);
744}
745
746proc lpGkDim(ideal G)
747"USAGE: lpGkDim(G); G an ideal in a letterplace ring
748RETURN: int
749PURPOSE: Determines the Gelfand Kirillov dimension of A/<G>
750@*       -1 means it is positive infinite
751ASSUME: - basering is a Letterplace ring
752@*      - G is a Groebner basis
753"
754{
755  G = lead(G);
756  G = simplify(G, 2+4+8);
757
758  // check special case 1
759  int l = 0;
760  for (int i = 1; i <= size(G); i++) {
761    // find the max degree in G
762    int d = deg(G[i]);
763    if (d > l) {
764      l = d;
765    }
766
767    // also if G is the whole ring return minus infinity
768    if (leadmonom(G[i]) == 1) {
769      ERROR("Gk-Dim not defined for 0-ring")
770    }
771    kill d;
772  } kill i;
773  // if longest word has length 1 we handle it as a special case
774  if (l == 1) {
775    int n = attrib(basering, "lV"); // variable count
776    int k = size(G);
777    if (k == n) { // V = {1} no edges
778      return(0);
779    }
780    if (k == n-1) { // V = {1} with loop
781      return(1);
782    }
783    if (k <= n-2) { // V = {1} with more than one loop
784      return(-1);
785    }
786  }
787
788  intmat UG = lpUfGraph(G);
789
790  // check special case 2
791  intmat zero[nrows(UG)][ncols(UG)];
792  if (UG == zero) {
793    return (0);
794  }
795
796  // check special case 3
797  UG = topologicalSort(UG);
798
799  if (imIsUpRightTriangle(UG)) {
800    UG = eliminateZerosUpTriangle(UG);
801    if (ncols(UG) == 0 || nrows(UG) == 0) { // when the diagonal was zero
802      return (0)
803    }
804    return(UfGraphURTNZDGrowth(UG));
805  }
806
807  // otherwise count cycles in the Ufnarovskij Graph
808  return(UfGraphGrowth(UG));
809}
810example
811{
812  "EXAMPLE:"; echo = 2;
813  ring r = 0,(x,y,z),dp;
814  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
815  R;
816  setring R; // sets basering to Letterplace ring
817  ideal I = z(1);//an example of infinite GK dimension
818  lpGkDim(I);
819  I = x(1),y(1),z(1); // gkDim = 0
820  lpGkDim(I);
821  I = x(1)*y(2), x(1)*z(2), z(1)*y(2), z(1)*z(2);//gkDim = 2
822  lpGkDim(I);
823}
824
825proc lpGlDimBound(ideal G)
826"USAGE: lpGlDimBound(I); I an ideal
827RETURN: int, an upper bound for the global dimension, -1 means infinity
828PURPOSE: computing an upper bound for the global dimension
829ASSUME: - basering is a Letterplace ring, G is a reduced Groebner Basis
830EXAMPLE: example lpGlDimBound; shows example
831NOTE: if I = LM(I), then the global dimension is equal the Gelfand
832    Kirillov dimension if it is finite
833    Global dimension should be 0 for A/G = K and 1 for A/G = K<x1...xn>
834"
835{
836  G = simplify(G,2); // remove zero generators
837  // NOTE: Gl should be 0 for A/G = K and 1 for A/G = K<x1...xn>
838  // G1 contains generators with single variable in LM
839  ideal G1;
840  for (int i = 1; i <= size(G); i++) {
841    if (ord(G[i]) < 2) { // single variable in LM
842      G1 = insertGenerator(G1,G[i]);
843    }
844  } kill i;
845  G1 = simplify(G1,2); // remove zero generators
846
847  // G = NF(G,G1)
848  for (int i = 1; i <= ncols(G); i++) { // do not use size() here
849    G[i] = lpNF(G[i],G1);
850  } kill i;
851  G = simplify(G,2); // remove zero generators
852
853  // delete variables in LM(G1) from the ring
854  def save = basering;
855  ring R = basering;
856  if (size(G1) > 0) {
857    while (size(G1) > 0) {
858      if (attrib(R, "lV") > 1) {
859        ring R = lpDelVar(lp2iv(G1[1])[1]);
860        ideal G1 = imap(save,G1);
861        G1 = simplify(G1, 2); // remove zero generators
862      } else {
863        // only the field is left (no variables)
864        return(0);
865      }
866    }
867    ideal G = imap(save, G); // put this here, because when save == R this call would make G = 0
868  }
869
870  // Li p. 184 if G = LM(G), then I = LM(I) and thus glDim = gkDim if it's finite
871  for (int i = 1; i <= size(G); i++) {
872    if (G[i] != lead(G[i])) {
873      break;
874    } else {
875      if (i == size(G)) { // if last iteration
876        int gkDim = lpGkDim(G);
877        if (gkDim >= 0) {
878          return (gkDim);
879        }
880        kill gkDim;
881      }
882    }
883  } kill i;
884
885  intmat GNC = lpGraphOfNChains(G);
886
887  // assuming GNC is connected
888
889  // TODO: maybe loop+cycle checking could be done more efficiently?
890  if (!imHasLoops(GNC) && imIsUpRightTriangle(topologicalSort(GNC))) {
891    // GNC is a DAG
892    intmat GNCk = GNC;
893    intmat zero[1][ncols(GNCk)];
894    int k = 1;
895    // while first row isn't empty
896    while (GNCk[1,1..(ncols(GNCk))] != zero[1,1..(ncols(zero))]) {
897      GNCk = GNCk * GNC;
898      k++;
899    }
900    // k-1 = number of edges in longest path starting from 1
901    return (k-1);
902  } else {
903    // GNC contains loops/cycles => there is always an n-chain
904    return (-1); // infinity
905  }
906}
907example
908{
909  "EXAMPLE:"; echo = 2;
910  ring r = 0,(x,y),dp;
911  def R = makeLetterplaceRing(5); // constructs a Letterplace ring
912  setring R; // sets basering to Letterplace ring
913  ideal G = x(1)*x(2), y(1)*y(2),x(1)*y(2)*x(3); // ideal G contains a
914  //Groebner basis
915  lpGlDimBound(G); // invokes procedure with Groebner basis G
916}
917
918static proc imHasLoops(intmat A) {
919  int n = ncols(A);
920  for (int i = 1; i < n; i++) {
921    if (A[i,i] == 1) {
922      return (1);
923    }
924  } kill i;
925  return (0);
926}
927
928static proc lpGraphOfNChains(ideal G) // G must be reduced
929{
930  list LG = lpId2ivLi(lead(G));
931  int n = attrib(basering, "lV");
932  int degbound = attrib(basering, "uptodeg");
933
934  list V;
935  for (int i = 0; i <= n; i++) {
936    V[i+1] = i; // add 1 and all variables
937  } kill i;
938  for (int i = 1; i <= size(LG); i++) {
939    intvec u = LG[i];
940    for (int j = 2; j <= size(u); j++) {
941      intvec v = u[j..size(u)];
942      if (!contains(V, v)) {
943        V = insert(V, v, size(V)); // add subword j..size
944      }
945      kill v;
946    } kill j;
947    kill u;
948  } kill i;
949  int nV = size(V);
950  intmat GNC[nV][nV]; // graph of n-chains
951
952  // for vertex 1
953  for (int i = 2; i <= n + 1; i++) {
954    GNC[1,i] = 1; // 1 has an edge to all variables
955  } kill i;
956
957  // for the other vertices
958  for (int i = 2; i <= nV; i++) {
959    for (int j = 2; j <= nV; j++) {
960      intvec uv = V[i],V[j];
961
962      if (contains(LG, uv)) {
963        GNC[i,j] = 1;
964      } else {
965        // Li p. 177
966        // search for a right divisor 'w' of uv in G
967        // then check if G doesn't divide the subword uv-1
968
969        // look for a right divisor in LG
970        for (int k = 1; k <= size(LG); k++) {
971          if (isSF(LG[k], uv)) {
972            // w = LG[k]
973            if(!ivdivides(LG, uv[1..(size(uv)-1)])) {
974              // G doesn't divide uv-1
975              GNC[i,j] = 1;
976              break;
977            }
978          }
979        } kill k;
980      }
981      kill uv;
982    } kill j;
983  } kill i;
984
985  return(GNC);
986}
987
988static proc contains(list L, def item)
989{
990  for (int i = 1; i <= size(L); i++) {
991    if (L[i] == item) {
992      return (1);
993    }
994  } kill i;
995  return (0);
996}
997
998/*proc lpSubstituteExpandRing(poly f, list s1, list s2) {*/
999/*int minDegBound = lpCalcSubstDegBound(f,s1,s2);*/
1000/**/
1001/*def R = basering; // curr lp ring*/
1002/*setring ORIGINALRING; // non lp ring TODO*/
1003/*def R1 = makeLetterplaceRing(minDegBound);*/
1004/*setring R1;*/
1005/**/
1006/*poly g = lpSubstitute(imap(R,f), imap(R,s1), imap(R,s2));*/
1007/**/
1008/*return (R1); // return the new ring*/
1009/*}*/
1010
1011proc lpSubstitute(poly f, ideal s1, ideal s2, list #)
1012"USAGE: lpSubstitute(f,s1,s2[,G]); f letterplace polynomial, s1 list (ideal) of variables
1013  to replace, s2 list (ideal) of polynomials to replace with, G optional ideal to
1014  reduce with.
1015RETURN: poly, the substituted polynomial
1016ASSUME: - basering is a Letterplace ring
1017      - G is a groebner basis,
1018      - the current ring has a sufficient degbound (can be calculated with
1019      lpCalcSubstDegBound())
1020EXAMPLE: example lpSubstitute; shows examples
1021"
1022{
1023  ideal G;
1024  if (size(#) > 0) {
1025    if (typeof(#[1])=="ideal") {
1026      G = #[1];
1027    }
1028  }
1029
1030  poly fs;
1031  for (int i = 1; i <= size(f); i++) {
1032    poly fis = leadcoef(f[i]);
1033    intvec ivfi = lp2iv(f[i]);
1034    for (int j = 1; j <= size(ivfi); j++) {
1035      int varindex = ivfi[j];
1036      int subindex = lpIndexOf(s1, var(varindex));
1037      if (subindex > 0) {
1038        s2[subindex] = lpNF(s2[subindex],G);
1039        fis = lpMult(fis, s2[subindex]);
1040      } else {
1041        fis = lpMult(fis, lpNF(iv2lp(varindex),G));
1042      }
1043      /*fis = lpNF(fis,G);*/
1044      kill varindex; kill subindex;
1045    } kill j;
1046    kill fis; kill ivfi;
1047    fs = fs + fis;
1048  } kill i;
1049  fs = lpNF(fs, G);
1050  return (fs);
1051}
1052example {
1053  //////// EXAMPLE A ////////
1054  LIB "fpadim.lib";
1055  ring r = 0,(x,y,z),dp;
1056  def R = makeLetterplaceRing(4);
1057  setring R;
1058
1059  ideal G = x(1)*y(2); // optional
1060
1061  poly f = 3*x(1)*x(2)+y(1)*x(2);
1062  ideal s1 = x(1), y(1);
1063  ideal s2 = y(1)*z(2)*z(3), x(1);
1064
1065  // the substitution probably needs a higher degbound
1066  int minDegBound = lpCalcSubstDegBounds(f,s1,s2);
1067  setring r;
1068  def R1 = makeLetterplaceRing(minDegBound);
1069  setring R1;
1070
1071  // the last parameter is optional
1072  lpSubstitute(imap(R,f), imap(R,s1), imap(R,s2), imap(R,G));
1073
1074  //////// EXAMPLE B ////////
1075  LIB "fpadim.lib";
1076  ring r = 0,(x,y,z),dp;
1077  def R = makeLetterplaceRing(4);
1078  setring R;
1079
1080  poly f = 3*x(1)*x(2)+y(1)*x(2);
1081  poly g = z(1)*x(2)+y(1);
1082  poly h = 7*x(1)*z(2)+x(1);
1083  ideal I = f,g,h;
1084  ideal s1 = x(1), y(1);
1085  ideal s2 = y(1)*z(2)*z(3), x(1);
1086
1087  int minDegBound = lpCalcSubstDegBounds(I,s1,s2);
1088  setring r;
1089  def R1 = makeLetterplaceRing(minDegBound);
1090  setring R1;
1091
1092  ideal I = imap(R,I);
1093  ideal s1 = imap(R,s1);
1094  ideal s2 = imap(R,s2);
1095  for (int i = 1; i <= size(I); i++) {
1096    lpSubstitute(I[i], s1, s2);
1097  }
1098}
1099
1100static proc lpIndexOf(ideal I, poly p) {
1101  for (int i = 1; i <= size(I); i++) {
1102    if (I[i] == p) {
1103      return (i);
1104    }
1105  } kill i;
1106  return (-1);
1107}
1108
1109static proc ivIndexOf(list L, intvec iv) {
1110  for (int i = 1; i <= size(L); i++) {
1111    if (L[i] == iv) {
1112      return (i);
1113    }
1114  } kill i;
1115  return (-1);
1116}
1117
1118
1119proc lpCalcSubstDegBound(poly f, ideal s1, ideal s2)
1120"USAGE: lpCalcSubstDegBound(f,s1,s2); f letterplace polynomial, s1 list (ideal) of variables
1121  to replace, s2 list (ideal) of polynomials to replace with
1122RETURN: int, the min degbound required to perform the substitution
1123ASSUME: - basering is a Letterplace ring
1124EXAMPLE: example lpSubstitute; shows examples
1125"
1126{
1127  int maxDegBound = 0;
1128  for (int i = 1; i <= size(f); i++) {
1129    intvec ivfi = lp2iv(f[i]);
1130    int tmpDegBound;
1131    for (int j = 1; j <= size(ivfi); j++) {
1132      int varindex = ivfi[j];
1133      int subindex = lpIndexOf(s1, var(varindex));
1134      if (subindex > 0) {
1135        tmpDegBound = tmpDegBound + deg(s2[subindex]);
1136      } else {
1137        tmpDegBound = tmpDegBound + 1;
1138      }
1139      kill varindex; kill subindex;
1140    } kill j;
1141    if (tmpDegBound > maxDegBound) {
1142      maxDegBound = tmpDegBound;
1143    }
1144    kill ivfi; kill tmpDegBound;
1145  } kill i;
1146
1147  // increase degbound by 50% when ideal is provided
1148  // needed for lpNF
1149  maxDegBound = maxDegBound + maxDegBound/2;
1150
1151  return (maxDegBound);
1152}
1153example {
1154  // see lpSubstitute()
1155}
1156
1157proc lpCalcSubstDegBounds(ideal I, ideal s1, ideal s2)
1158"USAGE: lpCalcSubstDegBounds(I,s1,s2); I list (ideal) of letterplace polynomials, s1 list (ideal)
1159  of variables to replace, s2 list (ideal) of polynomials to replace with
1160RETURN: int, the min degbound required to perform all of the substitutions
1161ASSUME: - basering is a Letterplace ring
1162EXAMPLE: example lpSubstitute; shows examples
1163NOTE: convenience method
1164"
1165{
1166  int maxDegBound = 0;
1167  for (int i = 1; i <= size(I); i++) {
1168    int tmpDegBound = lpCalcSubstDegBound(I[i], s1, s2, #);
1169    if (tmpDegBound > maxDegBound) {
1170      maxDegBound = tmpDegBound;
1171    }
1172    kill tmpDegBound;
1173  } kill i;
1174  return (maxDegBound);
1175}
1176example {
1177  // see lpSubstitute()
1178}
1179
1180static proc isSF(intvec S, intvec I)
1181"
1182PURPOSE:
1183checks, if a word S is a suffix of another word I
1184"
1185{
1186  int n = size(S);
1187  if (n <= 0 || S == 0) {return(1);}
1188  int m = size(I);
1189  if (m < n) {return(0);}
1190  intvec IS = I[(m-n+1)..m];
1191  if (IS == S) {return(1);}
1192  else {return(0);}
1193}
1194
1195static proc isIF(intvec IF, intvec I)
1196"
1197PURPOSE:
1198checks, if a word IF is an infix of another word I
1199"
1200{
1201  int n = size(IF);
1202  int m = size(I);
1203
1204  if (n <= 0 || IF == 0) {return(1);}
1205  if (m < n) {return(0);}
1206
1207  for (int i = 0; (n + i) <= m; i++){
1208    intvec IIF = I[(1 + i)..(n + i)];
1209    if (IIF == IF) {
1210      return(1);
1211    }
1212    kill IIF;
1213  } kill i;
1214  return(0);
1215}
1216
1217// removes a variable from a letterplace ring (a bit of a hack)
1218static proc lpDelVar(int index) {
1219  int lV = attrib(basering, "lV"); // number of variables in the main block
1220  int d = attrib(basering, "uptodeg"); // degree bround
1221  list LR = ringlist(basering);
1222
1223  if (!(index >= 1 && index <= lV)) { return (basering); } // invalid index
1224
1225  // remove frome the variable list
1226  for (int i = (d-1)*lV + index; i >= 1; i = i - lV) {
1227    LR[2] = delete(LR[2], i);
1228  } kill i;
1229
1230  // remove from a ordering
1231  intvec aiv = LR[3][1][2];
1232  aiv = aiv[1..(d*lV-d)];
1233  LR[3][1][2] = aiv;
1234
1235  // remove block orderings
1236  int del = (lV - index);
1237  int cnt = -1;
1238  for (int i = size(LR[3]); i >= 2; i--) {
1239    if (LR[3][i][2] != 0) {
1240      for (int j = size(LR[3][i][2]); j >= 1; j--) {
1241        cnt++; // next 1
1242        if (cnt%lV == del) {
1243          // delete
1244          if (size(LR[3][i][2]) > 1) { // if we have more than one element left, delete one
1245            LR[3][i][2] = delete(LR[3][i][2],j);
1246          } else { // otherwise delete the whole block
1247            LR[3] = delete(LR[3], i);
1248            break;
1249          }
1250        }
1251      } kill j;
1252    }
1253  } kill i;
1254
1255  def R = setLetterplaceAttributes(ring(LR),d,lV-1);
1256  return (R);
1257}
1258example
1259{
1260  "EXAMPLE:"; echo = 2;
1261  ring r = 0,(x,y,z),dp;
1262  def A = makeLetterplaceRing(3);
1263  setring A; A;
1264  def R = lpDelVar(2); setring R; R;
1265}
Note: See TracBrowser for help on using the repository browser.