source: git/Singular/LIB/fpaprops.lib @ 8a57d4

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