source: git/Singular/LIB/fpaprops.lib @ 6107b2

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