# source:git/Singular/LIB/freegb.lib@0610f0e

spielwiese
Last change on this file since 0610f0e was 0610f0e, checked in by Hans Schoenemann <hannes@…>, 13 years ago
• Property mode set to `100644`
File size: 67.2 KB
Line
1//////////////////////////////////////////////////////////////////////////////
2version="\$Id\$";
3category="Noncommutative";
4info="
5LIBRARY: freegb.lib   Compute two-sided Groebner bases in free algebras via letterplace
6AUTHOR: Viktor Levandovskyy,     levandov@math.rwth-aachen.de
7
8THEORY: See chapter 'LETTERPLACE' in the @sc{Singular} Manual.
9
10PROCEDURES:
11makeLetterplaceRing(d);    creates a ring with d blocks of shifted original variables
12freeGBasis(L, n);   compute two-sided Groebner basis of ideal, encoded via L, up to degree n
13setLetterplaceAttributes(R,d,b);  supplies ring R with the letterplace structure
14
15AUXILIARY PROCEDURES:
16
17lpMult(f,g);        letterplace multiplication of letterplace polynomials
18shiftPoly(p,i);    compute the i-th shift of letterplace polynomial p
19lpPower(f,n);     natural power of a letterplace polynomial
20lp2lstr(K, s);      convert letter-place ideal to a list of modules
21lst2str(L[, n]);     convert a list (of modules) into polynomials in free algebra
22mod2str(M[, n]); convert a module into a polynomial in free algebra
23vct2str(M[, n]);   convert a vector into a word in free algebra
24lieBracket(a,b[, N]);    compute Lie bracket ab-ba of two letterplace polynomials
25serreRelations(A,z);   compute the homogeneous part of Serre's relations associated to a generalized Cartan matrix A
26fullSerreRelations(A,N,C,P,d); compute the ideal of all Serre's relations associated to a generalized Cartan matrix A
27isVar(p);                   check whether p is a power of a single variable
28ademRelations(i,j);    compute the ideal of Adem relations for i<2j in char 0
29
31"
32
33// this library computes two-sided GB of an ideal
34// in a free associative algebra
35
36// a monomial is encoded via a vector V
37// where V[1] = coefficient
38// V[1+i] = the corresponding symbol
39
40LIB "qhmoduli.lib"; // for Max
41
42proc tstfreegb()
43{
44    /* tests all procs for consistency */
46  example makeLetterplaceRing;
47  example freeGBasis;
48  example setLetterplaceAttributes;
49  /* secondary */
50  example   lpMult;
51  example   shiftPoly;
52  example   lpPower;
53  example   lp2lstr;
54  example   lst2str;
55  example   mod2str;
56  example   vct2str;
57  example   lieBracket;
58  example   serreRelations;
59  example fullSerreRelations;
60  example   isVar;
62}
63
64proc setLetterplaceAttributes(def R, int uptodeg, int lV)
65"USAGE: setLetterplaceAttributes(R, d, b); R a ring, b,d integers
66RETURN: ring with special attributes set
67PURPOSE: sets attributes for a letterplace ring:
68@*      'isLetterplaceRing' = true, 'uptodeg' = d, 'lV' = b, where
69@*      'uptodeg' stands for the degree bound,
70@*      'lV' for the number of variables in the block 0.
71NOTE: Activate the resulting ring by using @code{setring}
72"
73{
74  if (uptodeg*lV != nvars(R))
75  {
76    ERROR("uptodeg and lV do not agree on the basering!");
77  }
78
79    // Set letterplace-specific attributes for the output ring!
80  attrib(R, "uptodeg", uptodeg);
81  attrib(R, "lV", lV);
82  attrib(R, "isLetterplaceRing", 1);
83  return (R);
84}
85example
86{
87  "EXAMPLE:"; echo = 2;
88  ring r = 0,(x(1),y(1),x(2),y(2),x(3),y(3),x(4),y(4)),dp;
89  def R = setLetterplaceAttributes(r, 4, 2); setring R;
90  attrib(R,"isLetterplaceRing");
91  lieBracket(x(1),y(1),2);
92}
93
94
95// obsolete?
96
97static proc lshift(module M, int s, string varing, def lpring)
98{
99  // FINALLY IMPLEMENTED AS A PART OT THE CODE
100  // shifts a polynomial from the ring R to s positions
101  // M lives in varing, the result in lpring
102  // to be run from varing
103  int i, j, k, sm, sv;
104  vector v;
105  //  execute("setring "+lpring);
106  setring lpring;
107  poly @@p;
108  ideal I;
109  execute("setring "+varing);
110  sm = ncols(M);
111  for (i=1; i<=s; i++)
112  {
113    // modules, e.g. free polynomials
114    for (j=1; j<=sm; j++)
115    {
116      //vectors, e.g. free monomials
117      v  = M[j];
118      sv = size(v);
119      sp = "@@p = @@p + ";
120      for (k=2; k<=sv; k++)
121      {
122        sp = sp + string(v[k])+"("+string(k-1+s)+")*";
123      }
124      sp = sp + string(v[1])+";"; // coef;
125      setring lpring;
126      //      execute("setring "+lpring);
127      execute(sp);
128      execute("setring "+varing);
129    }
130    setring lpring;
131    //    execute("setring "+lpring);
132    I = I,@@p;
133    @@p = 0;
134  }
135  setring lpring;
136  //execute("setring "+lpring);
137  export(I);
138  //  setring varing;
139  execute("setring "+varing);
140}
141
142static proc skip0(vector v)
143{
144  // skips zeros in a vector, producing another vector
145  if ( (v[1]==0) || (v==0) ) { return(vector(0)); }
146  int sv = nrows(v);
147  int sw = size(v);
148  if (sv == sw)
149  {
150    return(v);
151  }
152  int i;
153  int j=1;
154  vector w;
155  for (i=1; i<=sv; i++)
156  {
157    if (v[i] != 0)
158    {
159      w = w + v[i]*gen(j);
160      j++;
161    }
162  }
163  return(w);
164}
165
166proc lst2str(list L, list #)
167"USAGE:  lst2str(L[,n]);  L a list of modules, n an optional integer
168RETURN:  list (of strings)
169PURPOSE: convert a list (of modules) into polynomials in free algebra
170EXAMPLE: example lst2str; shows examples
171NOTE: if an optional integer is not 0, stars signs are used in multiplication
172"
173{
174  // returns a list of strings
175  // being sentences in words built from L
176  // if #[1] = 1, use * between generators
177  int useStar = 0;
178  if ( size(#)>0 )
179  {
180    if ( typeof(#[1]) != "int")
181    {
182      ERROR("Second argument of type int expected");
183    }
184    if (#[1])
185    {
186      useStar = 1;
187    }
188  }
189  int i;
190  int s    = size(L);
191  if (s<1) { return(list(""));}
192  list N;
193  for(i=1; i<=s; i++)
194  {
195    if ((typeof(L[i]) == "module") || (typeof(L[i]) == "matrix") )
196    {
197      N[i] = mod2str(L[i],useStar);
198    }
199    else
200    {
201      "module or matrix expected in the list";
202      return(N);
203    }
204  }
205  return(N);
206}
207example
208{
209  "EXAMPLE:"; echo = 2;
210  ring r = 0,(x,y,z),(dp(1),dp(2));
211  module M = [-1,x,y],[-7,y,y],[3,x,x];
212  module N = [1,x,y,x,y],[-2,y,x,y,x],[6,x,y,y,x,y];
213  list L; L[1] = M; L[2] = N;
214  lst2str(L);
215  lst2str(L[1],1);
216}
217
218
219proc mod2str(module M, list #)
220"USAGE:  mod2str(M[,n]);  M a module, n an optional integer
221RETURN:  string
222PURPOSE: convert a module into a polynomial in free algebra
223EXAMPLE: example mod2str; shows examples
224NOTE: if an optional integer is not 0, stars signs are used in multiplication
225"
226{
227  if (size(M)==0) { return(""); }
228  // returns a string
229  // a sentence in words built from M
230  // if #[1] = 1, use * between generators
231  int useStar = 0;
232  if ( size(#)>0 )
233  {
234    if ( typeof(#[1]) != "int")
235    {
236      ERROR("Second argument of type int expected");
237    }
238    if (#[1])
239    {
240      useStar = 1;
241    }
242  }
243  int i;
244  int s    = ncols(M);
245  string t;
246  string mp;
247  for(i=1; i<=s; i++)
248  {
249    mp = vct2str(M[i],useStar);
250    if (mp[1] == "-")
251    {
252      t = t + mp;
253    }
254    else
255    {
256      if (mp != "")
257      {
258         t = t + "+" + mp;
259      }
260    }
261  }
262  if (t[1]=="+")
263  {
264    t = t[2..size(t)]; // remove first "+"
265  }
266  return(t);
267}
268example
269{
270  "EXAMPLE:"; echo = 2;
271  ring r = 0,(x,y,z),(dp);
272  module M = [1,x,y,x,y],[-2,y,x,y,x],[6,x,y,y,x,y];
273  mod2str(M);
274  mod2str(M,1);
275}
276
277proc vct2str(vector v, list #)
278"USAGE:  vct2str(v[,n]);  v a vector, n an optional integer
279RETURN:  string
280PURPOSE: convert a vector into a word in free algebra
281EXAMPLE: example vct2str; shows examples
282NOTE: if an optional integer is not 0, stars signs are used in multiplication
283"
284{
285  if (v==0) { return(""); }
286  // if #[1] = 1, use * between generators
287  int useStar = 0;
288  if ( size(#)>0 )
289  {
290    if (#[1])
291    {
292      useStar = 1;
293    }
294  }
295  int ppl = printlevel-voice+2;
296  // for a word, encoded by v
297  // produces a string for it
298  v = skip0(v);
299  if (v==0) { return(string(""));}
301  int s = size(v);
302  string vs,vv,vp,err;
303  int i,j,p,q;
304  for (i=1; i<=s-1; i++)
305  {
306    p     = isVar(v[i+1]);
307    if (p==0)
308    {
309      err = "Error: monomial expected at nonzero position " + string(i+1);
310      ERROR(err+" in vct2str");
311      //      dbprint(ppl,err);
312      //      return("_");
313    }
314    if (p==1)
315    {
316      if (useStar && (size(vs) >0))       {   vs = vs + "*"; }
317        vs = vs + string(v[i+1]);
318    }
319    else //power
320    {
321      vv = string(v[i+1]);
322      q = find(vv,"^");
323      if (q==0)
324      {
325        q = find(vv,string(p));
326        if (q==0)
327        {
328          err = "error in find for string "+vv;
329          dbprint(ppl,err);
330          return("_");
331        }
332      }
333      // q>0
334      vp = vv[1..q-1];
335      for(j=1;j<=p;j++)
336      {
337         if (useStar && (size(vs) >0))       {   vs = vs + "*"; }
338         vs = vs + vp;
339      }
340    }
341  }
342  string scf;
343  if (cf == -1)
344  {
345    scf = "-";
346  }
347  else
348  {
349    scf = string(cf);
350    if ( (cf == 1) && (size(vs)>0) )
351    {
352      scf = "";
353    }
354  }
355  if (useStar && (size(scf) >0) && (scf!="-") )       {   scf = scf + "*"; }
356  vs = scf + vs;
357  return(vs);
358}
359example
360{
361  "EXAMPLE:"; echo = 2;
362  ring r = (0,a),(x,y3,z(1)),dp;
363  vector v = [-7,x,y3^4,x2,z(1)^3];
364  vct2str(v);
365  vct2str(v,1);
366  vector w = [-7a^5+6a,x,y3,y3,x,z(1),z(1)];
367  vct2str(w);
368  vct2str(w,1);
369}
370
371proc isVar(poly p)
372"USAGE:  isVar(p);  poly p
373RETURN:  int
374PURPOSE: check, whether leading monomial of p is a power of a single variable
375@* from the basering. Returns the exponent or 0 if p is multivariate.
376EXAMPLE: example isVar; shows examples
377"
378{
379  // checks whether p is a variable indeed
380  // if it's a power of a variable, returns the power
381  if (p==0) {  return(0); } //"p=0";
383  if ( (p-lead(p)) !=0 ) { return(0); } // "p-lm(p)>0";
385  int s = size(v);
386  int i=1;
387  int cnt = 0;
388  int pwr = 0;
389  for (i=1; i<=s; i++)
390  {
391    if (v[i] != 0)
392    {
393      cnt++;
394      pwr = v[i];
395    }
396  }
397  //  "cnt:";  cnt;
398  if (cnt==1) { return(pwr); }
399  else { return(0); }
400}
401example
402{
403  "EXAMPLE:"; echo = 2;
404  ring r = 0,(x,y),dp;
405  poly f = xy+1;
406  isVar(f);
407  poly g = y^3;
408  isVar(g);
409  poly h = 7*x^3;
410  isVar(h);
411  poly i = 1;
412  isVar(i);
413}
414
415// new conversion routines
416
417static proc id2words(ideal I, int d)
418{
419  // NOT FINISHED
420  // input: ideal I of polys in letter-place notation
421  // in the ring with d real vars
422  // output: the list of strings: associative words
423  // extract names of vars
424  int i,m,n; string s; string place = "(1)";
425  list lv;
426  for(i=1; i<=d; i++)
427  {
428    s = string(var(i));
429    // get rid of place
430    n = find(s, place);
431    if (n>0)
432    {
433      s = s[1..n-1];
434    }
435    lv[i] = s;
436  }
437  poly p,q;
438  for (i=1; i<=ncols(I); i++)
439  {
440    if (I[i] != 0)
441    {
442      p = I[i];
443      while (p!=0)
444      {
446      }
447    }
448  }
449
450  return(lv);
451}
452example
453{
454  "EXAMPLE:"; echo = 2;
455  ring r = 0,(x(1),y(1),z(1),x(2),y(2),z(2)),dp;
456  ideal I = x(1)*y(2) -z(1)*x(2);
457  id2words(I,3);
458}
459
460static proc mono2word(poly p, int d)
461{
462}
463
464// given the element -7xy^2x, it is represented as [-7,x,y^2,x] or as [-7,x,y,y,x]
465// use the orig ord on (x,y,z) and expand it blockwise to (x(i),y(i),z(i))
466
467// the correspondences:
468// monomial in K<x,y,z>    <<--->> vector in R
469// polynomial in K<x,y,z>  <<--->> list of vectors (matrix/module) in R
470// ideal in K<x,y,z>       <<--->> list of matrices/modules in R
471
472
473// 1. form a new ring
474// 2. NOP
475// 3. compute GB -> with the kernel stuff
476// 4. skip shifted elts (check that no such exist?)
477// 5. go back to orig vars, produce strings/modules
478// 6. return the result
479
480proc freeGBasis(list LM, int d)
481"USAGE:  freeGBasis(L, d);  L a list of modules, d an integer
482RETURN:  ring
483ASSUME: L has a special form. Namely, it is a list of modules, where
484@* each generator of every module stands for a monomial times coefficient in free algebra.
485@* In such a vector generator, the 1st entry is a nonzero coefficient from the ground field
486@* and each next entry hosts a variable from the basering.
487PURPOSE: compute the two-sided Groebner basis of an ideal, encoded by L
488@* in the free associative algebra, up to degree d
489NOTE: Apply @code{lst2str} to the output in order to obtain a better readable presentation
490EXAMPLE: example freeGBasis; shows examples
491"
492{
493  // d = up to degree, will be shifted to d+1
494  if (d<1) {"bad d"; return(0);}
495
496  int ppl = printlevel-voice+2;
497  string err = "";
498
499  int i,j,s;
500  def save = basering;
501  // determine max no of places in the input
502  int slm = size(LM); // numbers of polys in the ideal
503  int sm;
504  intvec iv;
505  module M;
506  for (i=1; i<=slm; i++)
507  {
508    // modules, e.g. free polynomials
509    M  = LM[i];
510    sm = ncols(M);
511    for (j=1; j<=sm; j++)
512    {
513      //vectors, e.g. free monomials
514      iv = iv, size(M[j])-1; // 1 place is reserved by the coeff
515    }
516  }
517  int D = Max(iv); // max size of input words
518  if (d<D) {"bad d"; return(LM);}
519  D = D + d-1;
520  //  D = d;
521  list LR  = ringlist(save);
522  list L, tmp;
523  L[1] = LR[1]; // ground field
524  L[4] = LR[4]; // quotient ideal
525  tmp  = LR[2]; // varnames
526  s = size(LR[2]);
527  for (i=1; i<=D; i++)
528  {
529    for (j=1; j<=s; j++)
530    {
531      tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")";
532    }
533  }
534  for (i=1; i<=s; i++)
535  {
536    tmp[i] = string(tmp[i])+"("+string(1)+")";
537  }
538  L[2] = tmp;
539  list OrigNames = LR[2];
540  // ordering: d blocks of the ord on r
541  // try to get whether the ord on r is blockord itself
542  s = size(LR[3]);
543  if (s==2)
544  {
545    // not a blockord, 1 block + module ord
546    tmp = LR[3][s]; // module ord
547    for (i=1; i<=D; i++)
548    {
549      LR[3][s-1+i] = LR[3][1];
550    }
551    LR[3][s+D] = tmp;
552  }
553  if (s>2)
554  {
555    // there are s-1 blocks
556    int nb = s-1;
557    tmp = LR[3][s]; // module ord
558    for (i=1; i<=D; i++)
559    {
560      for (j=1; j<=nb; j++)
561      {
562        LR[3][i*nb+j] = LR[3][j];
563      }
564    }
565    //    size(LR[3]);
566    LR[3][nb*(D+1)+1] = tmp;
567  }
568  L[3] = LR[3];
569  def @R = ring(L);
570  setring @R;
571  ideal I;
572  poly @p;
573  s = size(OrigNames);
574  //  "s:";s;
575  // convert LM to canonical vectors (no powers)
576  setring save;
577  kill M; // M was defined earlier
578  module M;
579  slm = size(LM); // numbers of polys in the ideal
580  int sv,k,l;
581  vector v;
582  //  poly p;
583  string sp;
584  setring @R;
585  poly @@p=0;
586  setring save;
587  for (l=1; l<=slm; l++)
588  {
589    // modules, e.g. free polynomials
590    M  = LM[l];
591    sm = ncols(M); // in intvec iv the sizes are stored
592    // modules, e.g. free polynomials
593    for (j=1; j<=sm; j++)
594    {
595      //vectors, e.g. free monomials
596      v  = M[j];
597      sv = size(v);
598      //        "sv:";sv;
599      sp = "@@p = @@p + ";
600      for (k=2; k<=sv; k++)
601      {
602        sp = sp + string(v[k])+"("+string(k-1)+")*";
603      }
604      sp = sp + string(v[1])+";"; // coef;
605      setring @R;
606      execute(sp);
607      setring save;
608    }
609    setring @R;
610    //      "@@p:"; @@p;
611    I = I,@@p;
612    @@p = 0;
613    setring save;
614  }
615  kill sp;
616  // 3. compute GB
617  setring @R;
618  dbprint(ppl,"computing GB");
619  ideal J = system("freegb",I,d,nvars(save));
620  //  ideal J = slimgb(I);
621  dbprint(ppl,J);
622  // 4. skip shifted elts
623  ideal K = select1(J,1..s); // s = size(OrigNames)
624  dbprint(ppl,K);
625  dbprint(ppl, "done with GB");
626  // K contains vars x(1),...z(1) = images of originals
627  // 5. go back to orig vars, produce strings/modules
628  if (K[1] == 0)
629  {
630    "no reasonable output, GB gives 0";
631    return(0);
632  }
633  int sk = size(K);
634  int sp, sx, a, b;
635  intvec x;
636  poly p,q;
637  poly pn;
638  // vars in 'save'
639  setring save;
640  module N;
641  list LN;
642  vector V;
643  poly pn;
644  // test and skip exponents >=2
645  setring @R;
646  for(i=1; i<=sk; i++)
647  {
648    p  = K[i];
649    while (p!=0)
650    {
652      //      "processing q:";q;
654      sx = size(x);
655      for(k=1; k<=sx; k++)
656      {
657        if ( x[k] >= 2 )
658        {
659          err = "skip: the value x[k] is " + string(x[k]);
660          dbprint(ppl,err);
661          //            return(0);
662          K[i] = 0;
663          p    = 0;
664          q    = 0;
665          break;
666        }
667      }
668      p  = p - q;
669    }
670  }
671  K  = simplify(K,2);
672  sk = size(K);
673  for(i=1; i<=sk; i++)
674  {
675    //    setring save;
676    //    V  = 0;
677    setring @R;
678    p  = K[i];
679    while (p!=0)
680    {
682      err =  "processing q:" + string(q);
683      dbprint(ppl,err);
685      sx = size(x);
687      setring save;
688      pn = imap(@R,pn);
689      V  = V + leadcoef(pn)*gen(1);
690      for(k=1; k<=sx; k++)
691      {
692        if (x[k] ==1)
693        {
694          a = k / s; // block number=a+1, a!=0
695          b = k % s; // remainder
696          //          printf("a: %s, b: %s",a,b);
697          if (b == 0)
698          {
699            // that is it's the last var in the block
700            b = s;
701            a = a-1;
702          }
703          V = V + var(b)*gen(a+2);
704        }
705//         else
706//         {
707//           printf("error: the value x[k] is %s", x[k]);
708//           return(0);
709//         }
710      }
711      err = "V: " + string(V);
712      dbprint(ppl,err);
713      //      printf("V: %s", string(V));
714      N = N,V;
715      V  = 0;
716      setring @R;
717      p  = p - q;
718      pn = 0;
719    }
720    setring save;
721    LN[i] = simplify(N,2);
722    N     = 0;
723  }
724  setring save;
725  return(LN);
726}
727example
728{
729  "EXAMPLE:"; echo = 2;
730  ring r = 0,(x,y,z),(dp(1),dp(2));
731  module M = [-1,x,y],[-7,y,y],[3,x,x]; // stands for free poly -xy - 7yy - 3xx
732  module N = [1,x,y,x],[-1,y,x,y]; // stands for free poly xyx - yxy
733  list L; L[1] = M; L[2] = N; // list of modules stands for an ideal in free algebra
734  lst2str(L); // list to string conversion of input polynomials
735  def U = freeGBasis(L,5); // 5 is the degree bound
736  lst2str(U);
737}
738
739static proc crs(list LM, int d)
740"USAGE:  crs(L, d);  L a list of modules, d an integer
741RETURN:  ring
742PURPOSE: create a ring and shift the ideal
743EXAMPLE: example crs; shows examples
744"
745{
746  // d = up to degree, will be shifted to d+1
747  if (d<1) {"bad d"; return(0);}
748
749  int ppl = printlevel-voice+2;
750  string err = "";
751
752  int i,j,s;
753  def save = basering;
754  // determine max no of places in the input
755  int slm = size(LM); // numbers of polys in the ideal
756  int sm;
757  intvec iv;
758  module M;
759  for (i=1; i<=slm; i++)
760  {
761    // modules, e.g. free polynomials
762    M  = LM[i];
763    sm = ncols(M);
764    for (j=1; j<=sm; j++)
765    {
766      //vectors, e.g. free monomials
767      iv = iv, size(M[j])-1; // 1 place is reserved by the coeff
768    }
769  }
770  int D = Max(iv); // max size of input words
771  if (d<D) {"bad d"; return(LM);}
772  D = D + d-1;
773  //  D = d;
774  list LR  = ringlist(save);
775  list L, tmp;
776  L[1] = LR[1]; // ground field
777  L[4] = LR[4]; // quotient ideal
778  tmp  = LR[2]; // varnames
779  s = size(LR[2]);
780  for (i=1; i<=D; i++)
781  {
782    for (j=1; j<=s; j++)
783    {
784      tmp[i*s+j] = string(tmp[j])+"("+string(i)+")";
785    }
786  }
787  for (i=1; i<=s; i++)
788  {
789    tmp[i] = string(tmp[i])+"("+string(0)+")";
790  }
791  L[2] = tmp;
792  list OrigNames = LR[2];
793  // ordering: d blocks of the ord on r
794  // try to get whether the ord on r is blockord itself
795  s = size(LR[3]);
796  if (s==2)
797  {
798    // not a blockord, 1 block + module ord
799    tmp = LR[3][s]; // module ord
800    for (i=1; i<=D; i++)
801    {
802      LR[3][s-1+i] = LR[3][1];
803    }
804    LR[3][s+D] = tmp;
805  }
806  if (s>2)
807  {
808    // there are s-1 blocks
809    int nb = s-1;
810    tmp = LR[3][s]; // module ord
811    for (i=1; i<=D; i++)
812    {
813      for (j=1; j<=nb; j++)
814      {
815        LR[3][i*nb+j] = LR[3][j];
816      }
817    }
818    //    size(LR[3]);
819    LR[3][nb*(D+1)+1] = tmp;
820  }
821  L[3] = LR[3];
822  def @R = ring(L);
823  setring @R;
824  ideal I;
825  poly @p;
826  s = size(OrigNames);
827  //  "s:";s;
828  // convert LM to canonical vectors (no powers)
829  setring save;
830  kill M; // M was defined earlier
831  module M;
832  slm = size(LM); // numbers of polys in the ideal
833  int sv,k,l;
834  vector v;
835  //  poly p;
836  string sp;
837  setring @R;
838  poly @@p=0;
839  setring save;
840  for (l=1; l<=slm; l++)
841  {
842    // modules, e.g. free polynomials
843    M  = LM[l];
844    sm = ncols(M); // in intvec iv the sizes are stored
845    for (i=0; i<=d-iv[l]; i++)
846    {
847      // modules, e.g. free polynomials
848      for (j=1; j<=sm; j++)
849      {
850        //vectors, e.g. free monomials
851        v  = M[j];
852        sv = size(v);
853        //        "sv:";sv;
854        sp = "@@p = @@p + ";
855        for (k=2; k<=sv; k++)
856        {
857          sp = sp + string(v[k])+"("+string(k-2+i)+")*";
858        }
859        sp = sp + string(v[1])+";"; // coef;
860        setring @R;
861        execute(sp);
862        setring save;
863      }
864      setring @R;
865      //      "@@p:"; @@p;
866      I = I,@@p;
867      @@p = 0;
868      setring save;
869    }
870  }
871  setring @R;
872  export I;
873  return(@R);
874}
875example
876{
877  "EXAMPLE:"; echo = 2;
878  ring r = 0,(x,y,z),(dp(1),dp(2));
879  module M = [-1,x,y],[-7,y,y],[3,x,x];
880  module N = [1,x,y,x],[-1,y,x,y];
881  list L; L[1] = M; L[2] = N;
882  lst2str(L);
883  def U = crs(L,5);
884  setring U; U;
885  I;
886}
887
888static proc polylen(ideal I)
889{
890  // returns the ideal of length of polys
891  int i;
892  intvec J;
893  number s = 0;
894  for(i=1;i<=ncols(I);i++)
895  {
896    J[i] = size(I[i]);
897    s = s + J[i];
898  }
899  printf("the sum of length %s",s);
900  //  print(s);
901  return(J);
902}
903
904// new: uniting both mLR1 (homog) and mLR2 (nonhomog)
905proc makeLetterplaceRing(int d, list #)
906"USAGE:  makeLetterplaceRing(d [,h]); d an integer, h an optional integer
907RETURN:  ring
908PURPOSE: creates a ring with the ordering, used in letterplace computations
909NOTE: if h is given an nonzero, the pure homogeneous letterplace block ordering will be used.
910EXAMPLE: example makeLetterplaceRing; shows examples
911"
912{
913  int use_old_mlr = 0;
914  if ( size(#)>0 )
915  {
916    if (( typeof(#[1]) == "int" ) || ( typeof(#[1]) == "poly" ) )
917    {
918      poly x = poly(#[1]);
919      if (x!=0)
920      {
921        use_old_mlr = 1;
922      }
923    }
924  }
925  if (use_old_mlr)
926  {
927    def @A = makeLetterplaceRing1(d);
928  }
929  else
930  {
931    def @A = makeLetterplaceRing2(d);
932  }
933  return(@A);
934}
935example
936{
937  "EXAMPLE:"; echo = 2;
938  ring r = 0,(x,y,z),(dp(1),dp(2));
939  def A = makeLetterplaceRing(2);
940  setring A;  A;
941  attrib(A,"isLetterplaceRing");
942  attrib(A,"uptodeg");  // degree bound
943  attrib(A,"lV"); // number of variables in the main block
944  setring r; def B = makeLetterplaceRing(2,1); // to compare:
945  setring B;  B;
946}
947
948//proc freegbRing(int d)
949proc makeLetterplaceRing1(int d)
950"USAGE:  makeLetterplaceRing1(d); d an integer
951RETURN:  ring
952PURPOSE: creates a ring with a special ordering, suitable for
953@* the use of homogeneous letterplace (d blocks of shifted original variables)
954EXAMPLE: example makeLetterplaceRing1; shows examples
955"
956{
957  // d = up to degree, will be shifted to d+1
958  if (d<1) {"bad d"; return(0);}
959
960  int uptodeg = d; int lV = nvars(basering);
961
962  int ppl = printlevel-voice+2;
963  string err = "";
964
965  int i,j,s;
966  def save = basering;
967  int D = d-1;
968  list LR  = ringlist(save);
969  list L, tmp;
970  L[1] = LR[1]; // ground field
971  L[4] = LR[4]; // quotient ideal
972  tmp  = LR[2]; // varnames
973  s = size(LR[2]);
974  for (i=1; i<=D; i++)
975  {
976    for (j=1; j<=s; j++)
977    {
978      tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")";
979    }
980  }
981  for (i=1; i<=s; i++)
982  {
983    tmp[i] = string(tmp[i])+"("+string(1)+")";
984  }
985  L[2] = tmp;
986  list OrigNames = LR[2];
987  // ordering: d blocks of the ord on r
988  // try to get whether the ord on r is blockord itself
989  // TODO: make L(2) ordering! exponent is maximally 2
990  s = size(LR[3]);
991  if (s==2)
992  {
993    // not a blockord, 1 block + module ord
994    tmp = LR[3][s]; // module ord
995    for (i=1; i<=D; i++)
996    {
997      LR[3][s-1+i] = LR[3][1];
998    }
999    LR[3][s+D] = tmp;
1000  }
1001  if (s>2)
1002  {
1003    // there are s-1 blocks
1004    int nb = s-1;
1005    tmp = LR[3][s]; // module ord
1006    for (i=1; i<=D; i++)
1007    {
1008      for (j=1; j<=nb; j++)
1009      {
1010        LR[3][i*nb+j] = LR[3][j];
1011      }
1012    }
1013    //    size(LR[3]);
1014    LR[3][nb*(D+1)+1] = tmp;
1015  }
1016  L[3] = LR[3];
1017  def @R = ring(L);
1018  //  setring @R;
1019  //  int uptodeg = d; int lV = nvars(basering); // were defined before
1020  def @@R = setLetterplaceAttributes(@R,uptodeg,lV);
1021  return (@@R);
1022}
1023example
1024{
1025  "EXAMPLE:"; echo = 2;
1026  ring r = 0,(x,y,z),(dp(1),dp(2));
1027  def A = makeLetterplaceRing1(2);
1028  setring A;
1029  A;
1030  attrib(A,"isLetterplaceRing");
1031  attrib(A,"uptodeg");  // degree bound
1032  attrib(A,"lV"); // number of variables in the main block
1033}
1034
1035proc makeLetterplaceRing2(int d)
1036"USAGE:  makeLetterplaceRing2(d); d an integer
1037RETURN:  ring
1038PURPOSE: creates a ring with a special ordering, suitable for
1039@* the use of non-homogeneous letterplace
1040NOTE: the matrix for the ordering looks as follows: first row is 1,1,...,1
1041@* then there come 'd' blocks of shifted original variables
1042EXAMPLE: example makeLetterplaceRing2; shows examples
1043"
1044{
1045  // d = up to degree, will be shifted to d+1
1046  if (d<1) {"bad d"; return(0);}
1047
1048  int uptodeg = d; int lV = nvars(basering);
1049
1050  int ppl = printlevel-voice+2;
1051  string err = "";
1052
1053  int i,j,s;
1054  def save = basering;
1055  int D = d-1;
1056  list LR  = ringlist(save);
1057  list L, tmp, tmp2, tmp3;
1058  L[1] = LR[1]; // ground field
1059  L[4] = LR[4]; // quotient ideal
1060  tmp  = LR[2]; // varnames
1061  s = size(LR[2]);
1062  for (i=1; i<=D; i++)
1063  {
1064    for (j=1; j<=s; j++)
1065    {
1066      tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")";
1067    }
1068  }
1069  for (i=1; i<=s; i++)
1070  {
1071    tmp[i] = string(tmp[i])+"("+string(1)+")";
1072  }
1073  L[2] = tmp;
1074  list OrigNames = LR[2];
1075  // ordering: one 1..1 a above
1076  // ordering: d blocks of the ord on r
1077  // try to get whether the ord on r is blockord itself
1078  // TODO: make L(2) ordering! exponent is maximally 2
1079  s = size(LR[3]);
1080  if (s==2)
1081  {
1082    // not a blockord, 1 block + module ord
1083    tmp = LR[3][s]; // module ord
1084    for (i=1; i<=d; i++)
1085    {
1086      LR[3][s-1+i] = LR[3][1];
1087    }
1088    //    LR[3][s+D] = tmp;
1089    LR[3][s+1+D] = tmp;
1090    LR[3][1] = list("a",intvec(1: int(d*lV))); // deg-ord
1091  }
1092  if (s>2)
1093  {
1094    // there are s-1 blocks
1095    int nb = s-1;
1096    tmp = LR[3][s]; // module ord to place at the very end
1097    tmp2 = LR[3]; tmp2 = tmp2[1..nb];
1098    tmp3[1] = list("a",intvec(1: int(d*lV))); // deg-ord, insert as the 1st
1099    for (i=1; i<=d; i++)
1100    {
1101      tmp3 = tmp3 + tmp2;
1102    }
1103    tmp3 = tmp3 + list(tmp);
1104    LR[3] = tmp3;
1105//     for (i=1; i<=d; i++)
1106//     {
1107//       for (j=1; j<=nb; j++)
1108//       {
1109//         //        LR[3][i*nb+j+1]= LR[3][j];
1110//         LR[3][i*nb+j+1]= tmp2[j];
1111//       }
1112//     }
1113//     //    size(LR[3]);
1114//     LR[3][(s-1)*d+2] = tmp;
1115//     LR[3] = list("a",intvec(1: int(d*lV))) + LR[3]; // deg-ord, insert as the 1st
1116    // remove everything behind nb*(D+1)+1 ?
1117    //    tmp = LR[3];
1118    //    LR[3] = tmp[1..size(tmp)-1];
1119  }
1120  L[3] = LR[3];
1121  def @R = ring(L);
1122  //  setring @R;
1123  //  int uptodeg = d; int lV = nvars(basering); // were defined before
1124  def @@R = setLetterplaceAttributes(@R,uptodeg,lV);
1125  return (@@R);
1126}
1127example
1128{
1129  "EXAMPLE:"; echo = 2;
1130  ring r = 0,(x,y,z),(dp(1),dp(2));
1131  def A = makeLetterplaceRing2(2);
1132  setring A;
1133  A;
1134  attrib(A,"isLetterplaceRing");
1135  attrib(A,"uptodeg");  // degree bound
1136  attrib(A,"lV"); // number of variables in the main block
1137}
1138
1139/* EXAMPLES:
1140
1141//static proc ex_shift()
1142{
1143  LIB "freegb.lib";
1144  ring r = 0,(x,y,z),(dp(1),dp(2));
1145  module M = [-1,x,y],[-7,y,y],[3,x,x];
1146  module N = [1,x,y,x],[-1,y,x,y];
1147  list L; L[1] = M; L[2] = N;
1148  lst2str(L);
1149  def U = crs(L,5);
1150  setring U; U;
1151  I;
1152  poly p = I[2]; // I[8];
1153  p;
1154  system("stest",p,7,7,3); // error -> the world is ok
1155  poly q1 = system("stest",p,1,7,3); //ok
1156  poly q6 = system("stest",p,6,7,3); //ok
1157  system("btest",p,3); //ok
1158  system("btest",q1,3); //ok
1159  system("btest",q6,3); //ok
1160}
1161
1162//static proc test_shrink()
1163{
1164  LIB "freegb.lib";
1165  ring r =0,(x,y,z),dp;
1166  int d = 5;
1167  def R = makeLetterplaceRing(d);
1168  setring R;
1169  poly p1 = x(1)*y(2)*z(3);
1170  poly p2 = x(1)*y(4)*z(5);
1171  poly p3 = x(1)*y(1)*z(3);
1172  poly p4 = x(1)*y(2)*z(2);
1173  poly p5 = x(3)*z(5);
1174  poly p6 = x(1)*y(1)*x(3)*z(5);
1175  poly p7 = x(1)*y(2)*x(3)*y(4)*z(5);
1176  poly p8 = p1+p2+p3+p4+p5 + p6 + p7;
1177  p1; system("shrinktest",p1,3);
1178  p2; system("shrinktest",p2,3);
1179  p3; system("shrinktest",p3,3);
1180  p4; system("shrinktest",p4,3);
1181  p5; system("shrinktest",p5,3);
1182  p6; system("shrinktest",p6,3);
1183  p7; system("shrinktest",p7,3);
1184  p8; system("shrinktest",p8,3);
1185  poly p9 = p1 + 2*p2 + 5*p5 + 7*p7;
1186  p9; system("shrinktest",p9,3);
1187}
1188
1189//static proc ex2()
1190{
1191  option(prot);
1192  LIB "freegb.lib";
1193  ring r = 0,(x,y),dp;
1194  module M = [-1,x,y],[3,x,x]; // 3x^2 - xy
1195  def U = freegb(M,7);
1196  lst2str(U);
1197}
1198
1199//static proc ex_nonhomog()
1200{
1201  option(prot);
1202  LIB "freegb.lib";
1203  ring r = 0,(x,y,h),dp;
1204  list L;
1205  module M;
1206  M = [-1,y,y],[1,x,x,x];  // x3-y2
1207  L[1] = M;
1208  M = [1,x,h],[-1,h,x];  // xh-hx
1209  L[2] = M;
1210  M = [1,y,h],[-1,h,y];  // yh-hy
1211  L[3] = M;
1212  def U = freegb(L,4);
1213  lst2str(U);
1214  // strange elements in the basis
1215}
1216
1217//static proc ex_nonhomog_comm()
1218{
1219  option(prot);
1220  LIB "freegb.lib";
1221  ring r = 0,(x,y),dp;
1222  module M = [-1,y,y],[1,x,x,x];
1223  def U = freegb(M,5);
1224  lst2str(U);
1225}
1226
1227//static proc ex_nonhomog_h()
1228{
1229  option(prot);
1230  LIB "freegb.lib";
1231  ring r = 0,(x,y,h),(a(1,1),dp);
1232  module M = [-1,y,y,h],[1,x,x,x]; // x3 - y2h
1233  def U = freegb(M,6);
1234  lst2str(U);
1235}
1236
1237//static proc ex_nonhomog_h2()
1238{
1239  option(prot);
1240  LIB "freegb.lib";
1241  ring r = 0,(x,y,h),(dp);
1242  list L;
1243  module M;
1244  M = [-1,y,y,h],[1,x,x,x]; // x3 - y2h
1245  L[1] = M;
1246  M = [1,x,h],[-1,h,x]; // xh - hx
1247  L[2] = M;
1248  M = [1,y,h],[-1,h,y]; // yh - hy
1249  L[3] = M;
1250  def U = freeGBasis(L,3);
1251  lst2str(U);
1253}
1254
1255
1256//static proc ex_nonhomog_3()
1257{
1258  option(prot);
1259  LIB "./freegb.lib";
1260  ring r = 0,(x,y,z),(dp);
1261  list L;
1262  module M;
1263  M = [1,z,y],[-1,x]; // zy - x
1264  L[1] = M;
1265  M = [1,z,x],[-1,y]; // zx - y
1266  L[2] = M;
1267  M = [1,y,x],[-1,z]; // yx - z
1268  L[3] = M;
1269  lst2str(L);
1270  list U = freegb(L,4);
1271  lst2str(U);
1273}
1274
1275//static proc ex_densep_2()
1276{
1277  option(prot);
1278  LIB "freegb.lib";
1279  ring r = (0,a,b,c),(x,y),(Dp); // deglex
1280  module M = [1,x,x], [a,x,y], [b,y,x], [c,y,y];
1281  lst2str(M);
1282  list U = freegb(M,5);
1283  lst2str(U);
1284  // a=b is important -> finite basis!!!
1285  module M = [1,x,x], [a,x,y], [a,y,x], [c,y,y];
1286  lst2str(M);
1287  list U = freegb(M,5);
1288  lst2str(U);
1289}
1290
1291// END COMMENTED EXAMPLES
1292
1293*/
1294
1295// 1. form a new ring
1296// 2. produce shifted generators
1297// 3. compute GB
1298// 4. skip shifted elts
1299// 5. go back to orig vars, produce strings/modules
1300// 6. return the result
1301
1302static proc freegbold(list LM, int d)
1303"USAGE:  freegbold(L, d);  L a list of modules, d an integer
1304RETURN:  ring
1305PURPOSE: compute the two-sided Groebner basis of an ideal, encoded by L in
1306the free associative algebra, up to degree d
1307EXAMPLE: example freegbold; shows examples
1308"
1309{
1310  // d = up to degree, will be shifted to d+1
1311  if (d<1) {"bad d"; return(0);}
1312
1313  int ppl = printlevel-voice+2;
1314  string err = "";
1315
1316  int i,j,s;
1317  def save = basering;
1318  // determine max no of places in the input
1319  int slm = size(LM); // numbers of polys in the ideal
1320  int sm;
1321  intvec iv;
1322  module M;
1323  for (i=1; i<=slm; i++)
1324  {
1325    // modules, e.g. free polynomials
1326    M  = LM[i];
1327    sm = ncols(M);
1328    for (j=1; j<=sm; j++)
1329    {
1330      //vectors, e.g. free monomials
1331      iv = iv, size(M[j])-1; // 1 place is reserved by the coeff
1332    }
1333  }
1334  int D = Max(iv); // max size of input words
1335  if (d<D) {"bad d"; return(LM);}
1336  D = D + d-1;
1337  //  D = d;
1338  list LR  = ringlist(save);
1339  list L, tmp;
1340  L[1] = LR[1]; // ground field
1341  L[4] = LR[4]; // quotient ideal
1342  tmp  = LR[2]; // varnames
1343  s = size(LR[2]);
1344  for (i=1; i<=D; i++)
1345  {
1346    for (j=1; j<=s; j++)
1347    {
1348      tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")";
1349    }
1350  }
1351  for (i=1; i<=s; i++)
1352  {
1353    tmp[i] = string(tmp[i])+"("+string(1)+")";
1354  }
1355  L[2] = tmp;
1356  list OrigNames = LR[2];
1357  // ordering: d blocks of the ord on r
1358  // try to get whether the ord on r is blockord itself
1359  // TODO: make L(2) ordering! exponent is maximally 2
1360  s = size(LR[3]);
1361  if (s==2)
1362  {
1363    // not a blockord, 1 block + module ord
1364    tmp = LR[3][s]; // module ord
1365    for (i=1; i<=D; i++)
1366    {
1367      LR[3][s-1+i] = LR[3][1];
1368    }
1369    LR[3][s+D] = tmp;
1370  }
1371  if (s>2)
1372  {
1373    // there are s-1 blocks
1374    int nb = s-1;
1375    tmp = LR[3][s]; // module ord
1376    for (i=1; i<=D; i++)
1377    {
1378      for (j=1; j<=nb; j++)
1379      {
1380        LR[3][i*nb+j] = LR[3][j];
1381      }
1382    }
1383    //    size(LR[3]);
1384    LR[3][nb*(D+1)+1] = tmp;
1385  }
1386  L[3] = LR[3];
1387  def @R = ring(L);
1388  setring @R;
1389  ideal I;
1390  poly @p;
1391  s = size(OrigNames);
1392  //  "s:";s;
1393  // convert LM to canonical vectors (no powers)
1394  setring save;
1395  kill M; // M was defined earlier
1396  module M;
1397  slm = size(LM); // numbers of polys in the ideal
1398  int sv,k,l;
1399  vector v;
1400  //  poly p;
1401  string sp;
1402  setring @R;
1403  poly @@p=0;
1404  setring save;
1405  for (l=1; l<=slm; l++)
1406  {
1407    // modules, e.g. free polynomials
1408    M  = LM[l];
1409    sm = ncols(M); // in intvec iv the sizes are stored
1410    for (i=0; i<=d-iv[l]; i++)
1411    {
1412      // modules, e.g. free polynomials
1413      for (j=1; j<=sm; j++)
1414      {
1415        //vectors, e.g. free monomials
1416        v  = M[j];
1417        sv = size(v);
1418        //        "sv:";sv;
1419        sp = "@@p = @@p + ";
1420        for (k=2; k<=sv; k++)
1421        {
1422          sp = sp + string(v[k])+"("+string(k-1+i)+")*";
1423        }
1424        sp = sp + string(v[1])+";"; // coef;
1425        setring @R;
1426        execute(sp);
1427        setring save;
1428      }
1429      setring @R;
1430      //      "@@p:"; @@p;
1431      I = I,@@p;
1432      @@p = 0;
1433      setring save;
1434    }
1435  }
1436  kill sp;
1437  // 3. compute GB
1438  setring @R;
1439  dbprint(ppl,"computing GB");
1440  //  ideal J = groebner(I);
1441  ideal J = slimgb(I);
1442  dbprint(ppl,J);
1443  // 4. skip shifted elts
1444  ideal K = select1(J,1..s); // s = size(OrigNames)
1445  dbprint(ppl,K);
1446  dbprint(ppl, "done with GB");
1447  // K contains vars x(1),...z(1) = images of originals
1448  // 5. go back to orig vars, produce strings/modules
1449  if (K[1] == 0)
1450  {
1451    "no reasonable output, GB gives 0";
1452    return(0);
1453  }
1454  int sk = size(K);
1455  int sp, sx, a, b;
1456  intvec x;
1457  poly p,q;
1458  poly pn;
1459  // vars in 'save'
1460  setring save;
1461  module N;
1462  list LN;
1463  vector V;
1464  poly pn;
1465  // test and skip exponents >=2
1466  setring @R;
1467  for(i=1; i<=sk; i++)
1468  {
1469    p  = K[i];
1470    while (p!=0)
1471    {
1473      //      "processing q:";q;
1475      sx = size(x);
1476      for(k=1; k<=sx; k++)
1477      {
1478        if ( x[k] >= 2 )
1479        {
1480          err = "skip: the value x[k] is " + string(x[k]);
1481          dbprint(ppl,err);
1482          //            return(0);
1483          K[i] = 0;
1484          p    = 0;
1485          q    = 0;
1486          break;
1487        }
1488      }
1489      p  = p - q;
1490    }
1491  }
1492  K  = simplify(K,2);
1493  sk = size(K);
1494  for(i=1; i<=sk; i++)
1495  {
1496    //    setring save;
1497    //    V  = 0;
1498    setring @R;
1499    p  = K[i];
1500    while (p!=0)
1501    {
1503      err =  "processing q:" + string(q);
1504      dbprint(ppl,err);
1506      sx = size(x);
1508      setring save;
1509      pn = imap(@R,pn);
1510      V  = V + leadcoef(pn)*gen(1);
1511      for(k=1; k<=sx; k++)
1512      {
1513        if (x[k] ==1)
1514        {
1515          a = k / s; // block number=a+1, a!=0
1516          b = k % s; // remainder
1517          //          printf("a: %s, b: %s",a,b);
1518          if (b == 0)
1519          {
1520            // that is it's the last var in the block
1521            b = s;
1522            a = a-1;
1523          }
1524          V = V + var(b)*gen(a+2);
1525        }
1526//         else
1527//         {
1528//           printf("error: the value x[k] is %s", x[k]);
1529//           return(0);
1530//         }
1531      }
1532      err = "V: " + string(V);
1533      dbprint(ppl,err);
1534      //      printf("V: %s", string(V));
1535      N = N,V;
1536      V  = 0;
1537      setring @R;
1538      p  = p - q;
1539      pn = 0;
1540    }
1541    setring save;
1542    LN[i] = simplify(N,2);
1543    N     = 0;
1544  }
1545  setring save;
1546  return(LN);
1547}
1548example
1549{
1550  "EXAMPLE:"; echo = 2;
1551  ring r = 0,(x,y,z),(dp(1),dp(2));
1552  module M = [-1,x,y],[-7,y,y],[3,x,x];
1553  module N = [1,x,y,x],[-1,y,x,y];
1554  list L; L[1] = M; L[2] = N;
1555  lst2str(L);
1556  def U = freegbold(L,5);
1557  lst2str(U);
1558}
1559
1560/* begin older procs and tests
1561
1562static proc sgb(ideal I, int d)
1563{
1564  // new code
1565  // map x_i to x_i(1) via map()
1566  //LIB "freegb.lib";
1567  def save = basering;
1568  //int d =7;// degree
1569  int nv = nvars(save);
1570  def R = makeLetterplaceRing(d);
1571  setring R;
1572  int i;
1573  ideal Imap;
1574  for (i=1; i<=nv; i++)
1575  {
1576    Imap[i] = var(i);
1577  }
1578  //ideal I = x(1)*y(2), y(1)*x(2)+z(1)*z(2);
1579  ideal I = x(1)*x(2),x(1)*y(2) + z(1)*x(2);
1580  option(prot);
1581  //option(teach);
1582  ideal J = system("freegb",I,d,nv);
1583}
1584
1585static proc checkCeq()
1586{
1587  ring r = 0,(x,y),Dp;
1588  def A = makeLetterplaceRing(4);
1589  setring A;
1590  A;
1591  // I = x2-xy
1592  ideal I = x(1)*x(2) - x(1)*y(2), x(2)*x(3) - x(2)*y(3), x(3)*x(4) - x(3)*y(4);
1593  ideal C = x(2)-x(1),x(3)-x(2),x(4)-x(3),y(2)-y(1),y(3)-y(2),y(4)-y(3);
1594  ideal K = I,C;
1595  groebner(K);
1596}
1597
1598static proc exHom1()
1599{
1601  // z*y - x, z*x - y, y*x - z
1602  LIB "freegb.lib";
1603  LIB "elim.lib";
1604  ring r = 0,(x,y,z,h),dp;
1605  list L;
1606  module M;
1607  M = [1,z,y],[-1,x,h]; // zy - xh
1608  L[1] = M;
1609  M = [1,z,x],[-1,y,h]; // zx - yh
1610  L[2] = M;
1611  M = [1,y,x],[-1,z,h]; // yx - zh
1612  L[3] = M;
1613  lst2str(L);
1614  def U = crs(L,4);
1615  setring U;
1616  I = I,
1617    y(2)*h(3)+z(2)*x(3),     y(3)*h(4)+z(3)*x(4),
1618    y(2)*x(3)-z(2)*h(3),     y(3)*x(4)-z(3)*h(4);
1619  I = simplify(I,2);
1620  ring r2 = 0,(x(0..4),y(0..4),z(0..4),h(0..4)),dp;
1621  ideal J = imap(U,I);
1622  //  ideal K = homog(J,h);
1623  option(redSB);
1624  option(redTail);
1625  ideal L = groebner(J); //(K);
1626  ideal LL = sat(L,ideal(h))[1];
1627  ideal M = subst(LL,h,1);
1628  M = simplify(M,2);
1629  setring U;
1630  ideal M = imap(r2,M);
1631  lst2str(U);
1632}
1633
1634static proc test1()
1635{
1636  LIB "freegb.lib";
1637  ring r = 0,(x,y),Dp;
1638  int d = 10; // degree
1639  def R = makeLetterplaceRing(d);
1640  setring R;
1641  ideal I = x(1)*x(2) - y(1)*y(2);
1642  option(prot);
1643  option(teach);
1644  ideal J = system("freegb",I,d,2);
1645  J;
1646}
1647
1648static proc test2()
1649{
1650  LIB "freegb.lib";
1651  ring r = 0,(x,y),Dp;
1652  int d = 10; // degree
1653  def R = makeLetterplaceRing(d);
1654  setring R;
1655  ideal I = x(1)*x(2) - x(1)*y(2);
1656  option(prot);
1657  option(teach);
1658  ideal J = system("freegb",I,d,2);
1659  J;
1660}
1661
1662static proc test3()
1663{
1664  LIB "freegb.lib";
1665  ring r = 0,(x,y,z),dp;
1666  int d =5; // degree
1667  def R = makeLetterplaceRing(d);
1668  setring R;
1669  ideal I = x(1)*y(2), y(1)*x(2)+z(1)*z(2);
1670  option(prot);
1671  option(teach);
1672  ideal J = system("freegb",I,d,3);
1673}
1674
1675static proc schur2-3()
1676{
1677  // nonhomog:
1678  //  h^4-10*h^2+9,f*e-e*f+h, h*2-e*h-2*e,h*f-f*h+2*f
1679  // homogenized with t
1680  //  h^4-10*h^2*t^2+9*t^4,f*e-e*f+h*t, h*2-e*h-2*e*t,h*f-f*h+2*f*t,
1681  // t*h - h*t, t*f - f*t, t*e - e*t
1682}
1683
1684end older procs and tests */
1685
1688RETURN:  ring (and exports ideal)
1689ASSUME: there are at least i+j variables in the basering
1690PURPOSE: compute the ideal of Adem relations for i<2j in characteristic 0
1691@*  the ideal is exported under the name AdemRel in the output ring
1693"
1694{
1695  // produces Adem relations for i<2j in char 0
1696  // assume: 0<i<2j
1697  // requires presence of vars up to i+j
1698  if ( (i<0) || (i >= 2*j) )
1699  {
1700    ERROR("arguments out of range"); return(0);
1701  }
1702  ring @r = 0,(s(i+j..0)),lp;
1703  poly p,q;
1704  number n;
1705  int ii = i div 2; int k;
1706  // k=0 => s(0)=1
1707  n = binomial(j-1,i);
1708  q = n*s(i+j)*s(0);
1709  //  printf("k=0, term=%s",q);
1710  p = p + q;
1711  for (k=1; k<= ii; k++)
1712  {
1713    n = binomial(j-k-1,i-2*k);
1714    q = n*s(i+j-k)*s(k);;
1715    //    printf("k=%s, term=%s",k,q);
1716    p = p + q;
1717  }
1720  return(@r);
1721}
1722example
1723{
1724  "EXAMPLE:"; echo = 2;
1726  setring A;
1728}
1729
1730/*
17311,1: 0
17321,2: s(3)*s(0) == s(3) -> def for s(3):=s(1)s(2)
17342,2: s(3)*s(1) == s(1)s(2)s(1)
17351,3: 0 ( since 2*s(4)*s(0) = 0 mod 2)
17372,3: s(5)*s(0)+s(4)*s(1) == s(5)+s(4)*s(1)
17383,2: 0
17393,3: s(5)*s(1)
17401,4: 3*s(5)*s(0) == s(5)  -> def for s(5):=s(1)*s(4)
17422,4: 3*s(6)*s(0)+s(5)*s(1) == s(6) + s(5)*s(1) == s(6) + s(1)*s(4)*s(1)
17444,3: s(5)*s(2)
17453,4: s(7)*s(0)+2*s(6)*s(1) == s(7) -> def for s(7):=s(3)*s(4)
17464,4: s(7)*s(1)+s(6)*s(2)
1747*/
1748
1749/* s1,s2:
1750s1*s1 =0, s2*s2 = s1*s2*s1
1751*/
1752
1753/*
1754try char 0:
1755s1,s2:
1756s1*s1 =0, s2*s2 = s1*s2*s1, s(1)*s(3)== s(1)*s(1)*s(3) == 0 = 2*s(4) ->def for s(4)
1757hence 2==0! only in char 2
1758 */
1759
1760  // Adem rels modulo 2 are interesting
1761
1762static proc stringpoly2lplace(string s)
1763{
1764  // decomposes sentence into terms
1765  s = replace(s,newline,""); // get rid of newlines
1766  s = replace(s," ",""); // get rid of empties
1767  //arith symbols: +,-
1768  // decompose into words with coeffs
1769  list LS;
1770  int i,j,ie,je,k,cnt;
1771  // s[1]="-" situation
1772  if (s[1]=="-")
1773  {
1774    LS = stringpoly2lplace(string(s[2..size(s)]));
1775    LS[1] = string("-"+string(LS[1]));
1776    return(LS);
1777  }
1778  i = find(s,"-",2);
1779  // i==1 might happen if the 1st symbol coeff is negative
1780  j = find(s,"+");
1781  list LL;
1782  if (i==j)
1783  {
1784    "return a monomial";
1785    // that is both are 0 -> s is a monomial
1786    LS[1] = s;
1787    return(LS);
1788  }
1789  if (i==0)
1790  {
1791    "i==0 situation";
1792    // no minuses at all => pluses only
1793    cnt++;
1794    LS[cnt] = string(s[1..j-1]);
1795    s = s[j+1..size(s)];
1796    while (s!= "")
1797    {
1798      j = find(s,"+");
1799      cnt++;
1800      if (j==0)
1801      {
1802        LS[cnt] = string(s);
1803        s = "";
1804      }
1805      else
1806      {
1807        LS[cnt] = string(s[1..j-1]);
1808        s = s[j+1..size(s)];
1809      }
1810    }
1811    return(LS);
1812  }
1813  if (j==0)
1814  {
1815    "j==0 situation";
1816    // no pluses at all except the lead coef => the rest are minuses only
1817    cnt++;
1818    LS[cnt] = string(s[1..i-1]);
1819    s = s[i..size(s)];
1820    while (s!= "")
1821    {
1822      i = find(s,"-",2);
1823      cnt++;
1824      if (i==0)
1825      {
1826        LS[cnt] = string(s);
1827        s = "";
1828      }
1829      else
1830      {
1831        LS[cnt] = string(s[1..i-1]);
1832        s = s[i..size(s)];
1833      }
1834    }
1835    return(LS);
1836  }
1837  // now i, j are nonzero
1838  if (i>j)
1839  {
1840    "i>j situation";
1841    // + comes first, at place j
1842    cnt++;
1843    //    "cnt:"; cnt; "j:"; j;
1844    LS[cnt] = string(s[1..j-1]);
1845    s = s[j+1..size(s)];
1846    LL = stringpoly2lplace(s);
1847    LS = LS + LL;
1848    kill LL;
1849    return(LS);
1850  }
1851  else
1852  {
1853    "j>i situation";
1854    // - might come first, at place i
1855    if (i>1)
1856    {
1857      cnt++;
1858      LS[cnt] = string(s[1..i-1]);
1859      s = s[i..size(s)];
1860    }
1861    else
1862    {
1863      // i==1->  minus at leadcoef
1864      ie = find(s,"-",i+1);
1865      je = find(s,"+",i+1);
1866      if (je == ie)
1867      {
1868         "ie=je situation";
1869        //monomial
1870        cnt++;
1871        LS[cnt] = s;
1872        return(LS);
1873      }
1874      if (je < ie)
1875      {
1876         "je<ie situation";
1877        // + comes first
1878        cnt++;
1879        LS[cnt] = s[1..je-1];
1880        s = s[je+1..size(s)];
1881      }
1882      else
1883      {
1884        // ie < je
1885         "ie<je situation";
1886        cnt++;
1887        LS[cnt] = s[1..ie-1];
1888        s = s[ie..size(s)];
1889      }
1890    }
1891    "going into recursion with "+s;
1892    LL = stringpoly2lplace(s);
1893    LS = LS + LL;
1894    return(LS);
1895  }
1896}
1897example
1898{
1899  "EXAMPLE:"; echo = 2;
1900  string s = "x*y+y*z+z*t"; // + only
1901  stringpoly2lplace(s);
1902  string s2 = "x*y - y*z-z*t*w*w"; // +1, - only
1903  stringpoly2lplace(s2);
1904  string s3 = "-x*y + y - 2*x +7*w*w*w";
1905  stringpoly2lplace(s3);
1906}
1907
1909{
1910  // adds places to the list of strings
1911  // according to their order in the list
1912  int s = size(L);
1913  int i;
1914  for (i=1; i<=s; i++)
1915  {
1916    if (typeof(L[i]) == "string")
1917    {
1918      L[i] = L[i] + "(" + string(i) + ")";
1919    }
1920    else
1921    {
1922      ERROR("entry of type string expected");
1923      return(0);
1924    }
1925  }
1926  return(L);
1927}
1928example
1929{
1930  "EXAMPLE:"; echo = 2;
1931  string a = "f1";   string b = "f2";
1932  list L = a,b,a;
1934}
1935
1936static proc sent2lplace(string s)
1937{
1938  // SENTence of words TO LetterPLACE
1939  list L =   stringpoly2lplace(s);
1940  int i; int ss = size(L);
1941  for(i=1; i<=ss; i++)
1942  {
1943    L[i] = str2lplace(L[i]);
1944  }
1945  return(L);
1946}
1947example
1948{
1949  "EXAMPLE:"; echo = 2;
1950  ring r = 0,(f2,f1),dp;
1951  string s = "f2*f1*f1 - 2*f1*f2*f1+ f1*f1*f2";
1952  sent2lplace(s);
1953}
1954
1955static proc testnumber(string s)
1956{
1957  string t;
1958  if (s[1]=="-")
1959  {
1960    // two situations: either there's a negative number
1961    t = s[2..size(s)];
1962    if (testnumber(t))
1963    {
1964      //a negative number
1965    }
1966    else
1967    {
1968      // a variable times -1
1969    }
1970    // or just a "-" for -1
1971  }
1972  t = "ring @r=(";
1973  t = t + charstr(basering)+"),";
1974  t = t + string(var(1))+",dp;";
1975  //  write(":w tstnum.tst",t);
1976  t = t+ "number @@Nn = " + s + ";"+"\$";
1977  write(":w tstnum.tst",t);
1978  string runsing = system("Singular");
1979  int k;
1980  t = runsing+ " -teq <tstnum.tst >tstnum.out";
1981  k = system("sh",t);
1982  if (k!=0)
1983  {
1984    ERROR("Problems running Singular");
1985  }
1986  int i = system("sh", "grep error tstnum.out > /dev/NULL");
1987  if (i!=0)
1988  {
1989    // no error: s is a number
1990    i = 1;
1991  }
1992  k = system("sh","rm tstnum.tst tstnum.out > /dev/NULL");
1993  return(i);
1994}
1995example
1996{
1997  "EXAMPLE:"; echo = 2;
1998  ring r = (0,a),x,dp;
1999  string s = "a^2+7*a-2";
2000  testnumber(s);
2001  s = "b+a";
2002  testnumber(s);
2003}
2004
2005static proc str2lplace(string s)
2006{
2007  // converts a word (monomial) with coeff into letter-place
2008  // string: coef*var1^exp1*var2^exp2*...varN^expN
2009  s = strpower2rep(s); // expand powers
2010  if (size(s)==0) { return(0); }
2011  int i,j,k,insC;
2012  string a,b,c,d,t;
2013  // 1. get coeff
2014  i = find(s,"*");
2015  if (i==0) { return(s); }
2016  list VN;
2017  c = s[1..i-1]; // incl. the case like (-a^2+1)
2018  int tn = testnumber(c);
2019  if (tn == 0)
2020  {
2021    // failed test
2022    if (c[1]=="-")
2023    {
2024      // two situations: either there's a negative number
2025      t = c[2..size(c)];
2026      if (testnumber(t))
2027      {
2028         //a negative number
2029        // nop here
2030      }
2031      else
2032      {
2033         // a variable times -1
2034          c = "-1";
2035          j++; VN[j] = t; //string(c[2..size(c)]);
2036          insC = 1;
2037      }
2038    }
2039    else
2040    {
2041      // just a variable with coeff 1
2042          j++; VN[j] = string(c);
2043          c = "1";
2044          insC = 1;
2045    }
2046  }
2047 // get vars
2048  t = s;
2049  //  t = s[i+1..size(s)];
2050  k = size(t); //j = 0;
2051  while (k>0)
2052  {
2053    t = t[i+1..size(t)]; //next part
2054    i = find(t,"*"); // next *
2055    if (i==0)
2056    {
2057      // last monomial
2058      j++;
2059      VN[j] = t;
2060      k = size(t);
2061      break;
2062    }
2063    b = t[1..i-1];
2064    //    print(b);
2065    j++;
2066    VN[j] = b;
2067    k = size(t);
2068  }
2070  VN[size(VN)+1] = string(c);
2071  return(VN);
2072}
2073example
2074{
2075  "EXAMPLE:"; echo = 2;
2076  ring r = (0,a),(f2,f1),dp;
2077  str2lplace("-2*f2^2*f1^2*f2");
2078  str2lplace("-f1*f2");
2079  str2lplace("(-a^2+7a)*f1*f2");
2080}
2081
2082static proc strpower2rep(string s)
2083{
2084  // makes x*x*x*x out of x^4 ., rep statys for repetitions
2085  // looks for "-" problem
2086  // exception: "-" as coeff
2087  string ex,t;
2088  int i,j,k;
2089
2090  i = find(s,"^"); // first ^
2091  if (i==0) { return(s); } // no ^ signs
2092
2093  if (s[1] == "-")
2094  {
2095    // either -coef or -1
2096    // got the coeff:
2097    i = find(s,"*");
2098    if (i==0)
2099    {
2100      // no *'s   => coef == -1 or s == -23
2101      i = size(s)+1;
2102    }
2103    t = string(s[2..i-1]); // without "-"
2104    if ( testnumber(t) )
2105    {
2106      // a good number
2107      t = strpower2rep(string(s[2..size(s)]));
2108      t = "-"+t;
2109      return(t);
2110    }
2111    else
2112    {
2113      // a variable
2114      t = strpower2rep(string(s[2..size(s)]));
2115      t = "-1*"+ t;
2116      return(t);
2117    }
2118  }
2119  // the case when leadcoef is a number in ()
2120  if (s[1] == "(")
2121  {
2122    i = find(s,")",2);    // must be nonzero
2123    t = s[2..i-1];
2124    if ( testnumber(t) )
2125    {
2126      // a good number
2127    }
2128    else {"strpower2rep: bad number as coef";}
2129    ex = string(s[i+2..size(s)]); // 2 because of *
2130    ex =  strpower2rep(ex);
2131    t = "("+t+")*"+ex;
2132    return(t);
2133  }
2134
2135  i = find(s,"^"); // first ^
2136  j = find(s,"*",i+1); // next * == end of ^
2137  if (j==0)
2138  {
2139    ex = s[i+1..size(s)];
2140  }
2141  else
2142  {
2143    ex = s[i+1..j-1];
2144  }
2145  execute("int @exp = " + ex + ";"); //@exp = exponent
2146  // got varname
2147  for (k=i-1; k>0; k--)
2148  {
2149    if (s[k] == "*") break;
2150  }
2151  string varn = s[k+1..i-1];
2152  //  "varn:";  varn;
2153  string pref;
2154  if (k>0)
2155  {
2156    pref = s[1..k]; // with * on the k-th place
2157  }
2158  //  "pref:";  pref;
2159  string suf;
2160  if ( (j>0) && (j+1 <= size(s)) )
2161  {
2162    suf = s[j+1..size(s)]; // without * on the 1st place
2163  }
2164  //  "suf:"; suf;
2165  string toins;
2166  for (k=1; k<=@exp; k++)
2167  {
2168    toins = toins + varn+"*";
2169  }
2170  //  "toins: ";  toins;
2171  if (size(suf) == 0)
2172  {
2173    toins = toins[1..size(toins)-1]; // get rid of trailing *
2174  }
2175  else
2176  {
2177    suf = strpower2rep(suf);
2178  }
2179  ex = pref + toins + suf;
2180  return(ex);
2181  //  return(strpower2rep(ex));
2182}
2183example
2184{
2185  "EXAMPLE:"; echo = 2;
2186  ring r = (0,a),(x,y,z,t),dp;
2187  strpower2rep("-x^4");
2188  strpower2rep("-2*x^4*y^3*z*t^2");
2189  strpower2rep("-a^2*x^4");
2190}
2191
2192proc lieBracket(poly a, poly b, list #)
2193"USAGE:  lieBracket(a,b[,N]); a,b letterplace polynomials, N an optional integer
2194RETURN:  poly
2195ASSUME: basering has a letterplace ring structure
2196PURPOSE: compute the Lie bracket [a,b] = ab - ba between letterplace polynomials
2197NOTE: if N>1 is specified, then the left normed bracket [a,[...[a,b]]]] is computed.
2198EXAMPLE: example lieBracket; shows examples
2199"
2200{
2201  if (lpAssumeViolation())
2202  {
2203    //    ERROR("Either 'uptodeg' or 'lV' global variables are not set!");
2204    ERROR("Incomplete Letterplace structure on the basering!");
2205  }
2206  // alias ppLiebr;
2207  //if int N is given compute [a,[...[a,b]]]] left normed bracket
2208  poly q;
2209  int N=1;
2210  if (size(#)>0)
2211  {
2212    if (typeof(#[1])=="int")
2213    {
2214      N = int(#[1]);
2215    }
2216  }
2217  if (N<=0) { return(q); }
2218  while (b!=0)
2219  {
2220    q = q + pmLiebr(a,lead(b));
2221    b = b - lead(b);
2222  }
2223  int i;
2224  if (N >1)
2225  {
2226    for(i=1; i<=N; i++)
2227    {
2228      q = lieBracket(a,q);
2229    }
2230  }
2231  return(q);
2232}
2233example
2234{
2235  "EXAMPLE:"; echo = 2;
2236  ring r = 0,(x(1),y(1),x(2),y(2),x(3),y(3),x(4),y(4)),dp;
2237  def R = setLetterplaceAttributes(r,4,2); // supply R with letterplace structure
2238  setring R;
2239  poly a = x(1)*y(2); poly b = y(1);
2240  lieBracket(a,b);
2241  lieBracket(x(1),y(1),2);
2242}
2243
2244static proc pmLiebr(poly a, poly b)
2245{
2246  //  a poly, b mono
2247  poly s;
2248  while (a!=0)
2249  {
2251    a = a - lead(a);
2252  }
2253  return(s);
2254}
2255
2256proc shiftPoly(poly a, int i)
2257"USAGE:  shiftPoly(p,i); p letterplace poly, i int
2258RETURN: poly
2259ASSUME: basering has letterplace ring structure
2260PURPOSE: compute the i-th shift of letterplace polynomial p
2261EXAMPLE: example shiftPoly; shows examples
2262"
2263{
2264  // shifts a monomial a by i
2265  // calls pLPshift(p,sh,uptodeg,lVblock);
2266  if (lpAssumeViolation())
2267  {
2268    ERROR("Incomplete Letterplace structure on the basering!");
2269  }
2270  int uptodeg = attrib(basering,"uptodeg");
2271  int lV = attrib(basering,"lV");
2272  if (deg(a) + i > uptodeg)
2273  {
2274    ERROR("degree bound violated by the shift!");
2275  }
2276  return(system("stest",a,i,uptodeg,lV));
2277}
2278example
2279{
2280  "EXAMPLE:"; echo = 2;
2281  ring r = 0,(x,y,z),dp;
2282  int uptodeg = 5; int lV = 3;
2283  def R = makeLetterplaceRing(uptodeg);
2284  setring R;
2285  poly f = x(1)*z(2)*y(3) - 2*z(1)*y(2) + 3*x(1);
2286  shiftPoly(f,1);
2287  shiftPoly(f,2);
2288}
2289
2290
2291static proc mmLiebr(poly a, poly b)
2292{
2293  // a,b, monomials
2296  int sa = deg(a);
2297  int sb = deg(b);
2298  poly v = a*shiftPoly(b,sa) - b*shiftPoly(a,sb);
2299  return(v);
2300}
2301
2302static proc test_shift()
2303{
2304  LIB "freegb.lib";
2305  ring r = 0,(a,b),dp;
2306  int d =5;
2307  def R = makeLetterplaceRing(d);
2308  setring R;
2309  int uptodeg = d;
2310  int lV = 2;
2311  def R = setLetterplaceAttributes(r,uptodeg,2); // supply R with letterplace structure
2312  setring R;
2313  poly p = mmLiebr(a(1),b(1));
2314  poly p = lieBracket(a(1),b(1));
2315}
2316
2317proc serreRelations(intmat A, int zu)
2318"USAGE:  serreRelations(A,z); A an intmat, z an int
2319RETURN:  ideal
2320ASSUME: basering has a letterplace ring structure and
2321@*          A is a generalized Cartan matrix with integer entries
2322PURPOSE: compute the ideal of Serre's relations associated to A
2323EXAMPLE: example serreRelations; shows examples
2324"
2325{
2326  // zu = 1 -> with commutators [f_i,f_j]; zu == 0 without them
2327  // suppose that A is cartan matrix
2328  // then Serre's relations are
2329  // (ad f_j)^{1-A_{ij}} ( f_i)
2330  int ppl = printlevel-voice+2;
2331  int n = ncols(A); // hence n variables
2332  int i,j,k,el;
2333  poly p,q;
2334  ideal I;
2335  for (i=1; i<=n; i++)
2336  {
2337    for (j=1; j<=n; j++)
2338    {
2339      el = 1 - A[i,j];
2340      //     printf("i:%s, j: %s, l: %s",i,j,l);
2341      dbprint(ppl,"i, j, l: ",i,j,el);
2342      //      if ((i!=j) && (l >0))
2343      //      if ( (i!=j) &&  ( ((zu ==0) &&  (l >=2)) || ((zu ==1) &&  (l >=1)) ) )
2344      if ((i!=j) && (el >0))
2345      {
2346        q = lieBracket(var(j),var(i));
2347        dbprint(ppl,"first bracket: ",q);
2348        //        if (l >=2)
2349        //        {
2350          for (k=1; k<=el-1; k++)
2351          {
2352            q = lieBracket(var(j),q);
2353            dbprint(ppl,"further bracket:",q);
2354          }
2355          //        }
2356      }
2357      if (q!=0) { I = I,q; q=0;}
2358    }
2359  }
2360  I = simplify(I,2);
2361  return(I);
2362}
2363example
2364{
2365  "EXAMPLE:"; echo = 2;
2366  intmat A[3][3] =
2367    2, -1, 0,
2368    -1, 2, -3,
2369    0, -1, 2; // G^1_2 Cartan matrix
2370  ring r = 0,(f1,f2,f3),dp;
2371  int uptodeg = 5;
2372  def R = makeLetterplaceRing(uptodeg);
2373  setring R;
2374  ideal I = serreRelations(A,1); I = simplify(I,1+2+8);
2375  I;
2376}
2377
2378/* setup for older example:
2379  intmat A[2][2] = 2, -1, -1, 2; // sl_3 == A_2
2380  ring r = 0,(f1,f2),dp;
2381  int uptodeg = 5; int lV = 2;
2382*/
2383
2384proc fullSerreRelations(intmat A, ideal rNegative, ideal rCartan, ideal rPositive, int uptodeg)
2385"USAGE:  fullSerreRelations(A,N,C,P,d); A an intmat, N,C,P ideals, d an int
2386RETURN:  ring (and ideal)
2387PURPOSE: compute the inhomogeneous Serre's relations associated to A in given variable names
2388ASSUME: three ideals in the input are of the same sizes and contain merely variables
2389@* which are interpreted as follows: N resp. P stand for negative resp. positive roots,
2390@* C stand for Cartan elements. d is the degree bound for letterplace ring, which will be returned.
2391@* The matrix A is a generalized Cartan matrix with integer entries
2392@* The result is the ideal called 'fsRel' in the returned ring.
2393EXAMPLE: example fullSerreRelations; shows examples
2394"
2395{
2396  /* SerreRels on rNeg and rPos plus Cartans etc. */
2397  int ppl = printlevel -voice+2;
2398  /* ideals must be written in variables: assume each term is of degree 1 */
2399  int i,j,k;
2400  int N = nvars(basering);
2401  def save = basering;
2402  int comFlag = 0;
2403  /* assume:  (size(rNegative) == size(rPositive)) */
2404  /* assume:  (size(rNegative) == size(rCartan)) i.e. nonsimple Cartans */
2405  if ( (size(rNegative) != size(rPositive)) || (size(rNegative) != size(rCartan)) )
2406  {
2407    ERROR("All input ideals must be of the same size");
2408  }
2409
2410//   if (size(rNegative) != size(rPositive))
2411//   {
2412//     ERROR("The 1st and the 3rd input ideals must be of the same size");
2413//   }
2414
2415  /* assume:  2*size(rNegative) + size(rCartan) >= nvars(basering) */
2416  i = 2*size(rNegative) + size(rCartan);
2417  if (i>N)
2418  {
2419    ERROR("The total number of elements in input ideals must not exceed the dimension of the ground ring");
2420  }
2421  if (i < N)
2422  {
2423    comFlag = N-i; // so many elements will commute
2424    "Warning: some elements will be treated as mutually commuting";
2425  }
2426  /* extract varnames from input ideals */
2427  intvec iNeg = varIdeal2intvec(rNegative);
2428  intvec iCartan = varIdeal2intvec(rCartan);
2429  intvec iPos = varIdeal2intvec(rPositive);
2430  /* for each vector in rNeg and rPositive, go into the corr. ring and create SerreRels */
2431  /* rNegative: */
2432  list L = ringlist(save);
2433  def LPsave = makeLetterplaceRing2(uptodeg); setring save;
2434  list LNEG = L; list tmp;
2435  /* L[1] field as is; L[2] vars: a subset; L[3] ordering: dp, L[4] as is */
2436  for (i=1; i<=size(iNeg); i++)
2437  {
2438    tmp[i] = string(var(iNeg[i]));
2439  }
2440  LNEG[2] = tmp; LNEG[3] = list(list("dp",intvec(1:size(iNeg))), list("C",0));
2441  def RNEG = ring(LNEG); setring RNEG;
2442  def RRNEG = makeLetterplaceRing2(uptodeg);
2443  setring RRNEG;
2444  ideal I = serreRelations(A,1); I = simplify(I,1+2+8);
2445  setring LPsave;
2446  ideal srNeg = imap(RRNEG,I);
2447  dbprint(ppl,"0-1 ideal of negative relations is ready");
2448  dbprint(ppl-1,srNeg);
2449  setring save; kill L,tmp,RRNEG,RNEG, LNEG;
2450  /* rPositive: */
2451  list L = ringlist(save);
2452  list LPOS = L; list tmp;
2453  /* L[1] field as is; L[2] vars: a subset; L[3] ordering: dp, L[4] as is */
2454  for (i=1; i<=size(iPos); i++)
2455  {
2456    tmp[i] = string(var(iPos[i]));
2457  }
2458  LPOS[2] = tmp; LPOS[3] = list(list("dp",intvec(1:size(iPos))), list("C",0));
2459  def RPOS = ring(LPOS); setring RPOS;
2460  def RRPOS = makeLetterplaceRing2(uptodeg);
2461  setring RRPOS;
2462  ideal I = serreRelations(A,1); I = simplify(I,1+2+8);
2463  setring LPsave;
2464  ideal srPos = imap(RRPOS,I);
2465  dbprint(ppl,"0-2 ideal of positive relations is ready");
2466  dbprint(ppl-1,srPos);
2467  setring save; kill L,tmp,RRPOS,RPOS, LPOS;
2468  string sMap = "ideal Mmap =";
2469  for (i=1; i<=nvars(save); i++)
2470  {
2471    sMap = sMap + string(var(i)) +"(1),";
2472  }
2473  sMap[size(sMap)] = ";";
2474  /* cartans: h_j h_i = h_i h_j */
2475  setring LPsave;
2476  ideal ComCartan;
2477  for (i=1; i<size(iCartan); i++)
2478  {
2479    for (j=i+1; j<=size(iCartan); j++)
2480    {
2481      ComCartan =  ComCartan + lieBracket(var(iCartan[j]),var(iCartan[i]));
2482    }
2483  }
2484  ComCartan = simplify(ComCartan,1+2+8);
2485  execute(sMap); // defines an ideal Mmap
2486  map F = save, Mmap;
2487  dbprint(ppl,"1. commuting Cartans: ");
2488  dbprint(ppl-1,ComCartan);
2489  /* [e_i, f_j] =0 if i<>j */
2490  ideal ComPosNeg; // assume: #Neg=#Pos
2491  for (i=1; i<size(iPos); i++)
2492  {
2493    for (j=1; j<=size(iPos); j++)
2494    {
2495      if (j !=i)
2496      {
2497        ComPosNeg =  ComPosNeg + lieBracket(var(iPos[i]),var(iNeg[j]));
2498        ComPosNeg =  ComPosNeg + lieBracket(var(iPos[j]),var(iNeg[i]));
2499      }
2500    }
2501  }
2502  ComPosNeg = simplify(ComPosNeg,1+2+8);
2503  dbprint(ppl,"2. commuting Positive and Negative:");
2504  dbprint(ppl-1,ComPosNeg);
2505  /* [e_i, f_i] = h_i */
2506  poly tempo;
2507  for (i=1; i<=size(iCartan); i++)
2508  {
2509    tempo = lieBracket(var(iPos[i]),var(iNeg[i])) - var(iCartan[i]);
2510    ComPosNeg =  ComPosNeg + tempo;
2511  }
2512  //  ComPosNeg = simplify(ComPosNeg,1+2+8);
2513  dbprint(ppl,"3. added sl2 triples [e_i,f_i]=h_i");
2514  dbprint(ppl-1,ComPosNeg);
2515
2516  /* [h_i, e_j] = A_ij e_j */
2517  /* [h_i, f_j] = -A_ij f_j */
2518  ideal ActCartan; // assume: #Neg=#Pos
2519  for (i=1; i<=size(iCartan); i++)
2520  {
2521    for (j=1; j<=size(iCartan); j++)
2522    {
2523      tempo = lieBracket(var(iCartan[i]),var(iPos[j])) - A[i,j]*var(iPos[j]);
2524      ActCartan = ActCartan + tempo;
2525      tempo = lieBracket(var(iCartan[i]),var(iNeg[j])) + A[i,j]*var(iNeg[j]);
2526      ActCartan = ActCartan + tempo;
2527    }
2528  }
2529  ActCartan = simplify(ActCartan,1+2+8);
2530  dbprint(ppl,"4. actions of Cartan:");
2531  dbprint(ppl-1, ActCartan);
2532
2533  /* final part: prepare the output */
2534  setring LPsave;
2535  ideal fsRel = srNeg, srPos, ComPosNeg, ComCartan, ActCartan;
2536  export fsRel;
2537  setring save;
2538  return(LPsave);
2539}
2540example
2541{
2542  "EXAMPLE:"; echo = 2;
2543  intmat A[2][2] =
2544    2, -1,
2545    -1, 2; // A_2 = sl_3 Cartan matrix
2546  ring r = 0,(f1,f2,h1,h2,e1,e2),dp;
2547  ideal negroots = f1,f2; ideal cartans = h1,h2; ideal posroots = e1,e2;
2548  int uptodeg = 5;
2549  def RS = fullSerreRelations(A,negroots,cartans,posroots,uptodeg);
2550  setring RS; fsRel;
2551}
2552
2553proc varIdeal2intvec(ideal I)
2554{
2555  /* assume1:  input ideal is a list of variables of the ground ring */
2556  int i,j; intvec V;
2557  for (i=1; i<= size(I); i++)
2558  {
2559    j = univariate(I[i]);
2560    if (j<=0)
2561    {
2562      ERROR("input ideal must contain only variables");
2563    }
2564    V[i] = j;
2565  }
2566  dbprint(printlevel-voice+2,V);
2567  /* now we make a smaller list of non-repeating entries */
2568  ideal iW = simplify(ideal(V),2+4); // no zeros, no repetitions
2569  if (size(iW) < size(V))
2570  {
2571    /* extract intvec from iW */
2572    intvec inW;
2573    for(j=1; j<=size(iW); j++)
2574    {
2576    }
2577    return(inW);
2578  }
2579  return(V);
2580}
2581example
2582{
2583  "EXAMPLE:"; echo = 2;
2584  ring r = 0,(x,y,z),dp;
2585  ideal I = x,z;
2586  varIdeal2intvec(I);
2587  varIdeal2intvec(ideal(x2,y^3,x+1));
2588  varIdeal2intvec(ideal(x*y,y,x+1));
2589}
2590
2591proc lp2lstr(ideal K, def save)
2592"USAGE:  lp2lstr(K,s); K an ideal, s a ring name
2593RETURN:  nothing (exports object @LN into the ring named s)
2594ASSUME: basering has a letterplace ring structure
2595PURPOSE: converts letterplace ideal to list of modules
2596NOTE: useful as preprocessing to 'lst2str'
2597EXAMPLE: example lp2lstr; shows examples
2598"
2599{
2600  def @R = basering;
2601  string err;
2602  int s = nvars(save);
2603  int i,j,k;
2604    // K contains vars x(1),...z(1) = images of originals
2605  // 5. go back to orig vars, produce strings/modules
2606  int sk = size(K);
2607  int sp, sx, a, b;
2608  intvec x;
2609  poly p,q;
2610  poly pn;
2611  // vars in 'save'
2612  setring save;
2613  module N;
2614  list LN;
2615  vector V;
2616  poly pn;
2617  // test and skip exponents >=2
2618  setring @R;
2619  for(i=1; i<=sk; i++)
2620  {
2621    p  = K[i];
2622    while (p!=0)
2623    {
2625      //      "processing q:";q;
2627      sx = size(x);
2628      for(k=1; k<=sx; k++)
2629      {
2630        if ( x[k] >= 2 )
2631        {
2632          err = "skip: the value x[k] is " + string(x[k]);
2633          dbprint(ppl,err);
2634          //            return(0);
2635          K[i] = 0;
2636          p    = 0;
2637          q    = 0;
2638          break;
2639        }
2640      }
2641      p  = p - q;
2642    }
2643  }
2644  K  = simplify(K,2);
2645  sk = size(K);
2646  for(i=1; i<=sk; i++)
2647  {
2648    //    setring save;
2649    //    V  = 0;
2650    setring @R;
2651    p  = K[i];
2652    while (p!=0)
2653    {
2655      err =  "processing q:" + string(q);
2656      dbprint(ppl,err);
2658      sx = size(x);
2660      setring save;
2661      pn = imap(@R,pn);
2662      V  = V + leadcoef(pn)*gen(1);
2663      for(k=1; k<=sx; k++)
2664      {
2665        if (x[k] ==1)
2666        {
2667          a = k / s; // block number=a+1, a!=0
2668          b = k % s; // remainder
2669          //          printf("a: %s, b: %s",a,b);
2670          if (b == 0)
2671          {
2672            // that is it's the last var in the block
2673            b = s;
2674            a = a-1;
2675          }
2676          V = V + var(b)*gen(a+2);
2677        }
2678      }
2679      err = "V: " + string(V);
2680      dbprint(ppl,err);
2681      //      printf("V: %s", string(V));
2682      N = N,V;
2683      V  = 0;
2684      setring @R;
2685      p  = p - q;
2686      pn = 0;
2687    }
2688    setring save;
2689    LN[i] = simplify(N,2);
2690    N     = 0;
2691  }
2692  setring save;
2693  list @LN = LN;
2694  export @LN;
2695  //  return(LN);
2696}
2697example
2698{
2699  "EXAMPLE:"; echo = 2;
2700  intmat A[2][2] = 2, -1, -1, 2; // sl_3 == A_2
2701  ring r = 0,(f1,f2),dp;
2702  def R = makeLetterplaceRing(3);
2703  setring R;
2704  ideal I = serreRelations(A,1);
2705  lp2lstr(I,r);
2706  setring r;
2707  lst2str(@LN,1);
2708}
2709
2710static proc strList2poly(list L)
2711{
2712  //  list L comes from sent2lplace (which takes a polynomial as input)
2713  // each entry of L is a sublist with the coef on the last place
2714  int s = size(L); int t;
2715  int i,j;
2716  list M;
2717  poly p,q;
2718  string Q;
2719  for(i=1; i<=s; i++)
2720  {
2721    M = L[i];
2722    t = size(M);
2723    //    q = M[t]; // a constant
2724    Q = string(M[t]);
2725    for(j=1; j<t; j++)
2726    {
2727      //      q = q*M[j];
2728      Q = Q+"*"+string(M[j]);
2729    }
2730    execute("q="+Q+";");
2731    //    q;
2732    p = p + q;
2733  }
2734  kill Q;
2735  return(p);
2736}
2737example
2738{
2739  "EXAMPLE:"; echo = 2;
2740  ring r =0,(x,y,z,t),Dp;
2741  def A = makeLetterplaceRing(4);
2742  setring A;
2743  string t = "-2*y*z*y*z + y*t*z*z - z*x*x*y  + 2*z*y*z*y";
2744  list L = sent2lplace(t);
2745  L;
2746  poly p = strList2poly(L);
2747  p;
2748}
2749
2750static proc file2lplace(string fname)
2751"USAGE:  file2lplace(fnm);  fnm a string
2752RETURN:  ideal
2753PURPOSE: convert the contents of the file fnm into ideal of polynomials in free algebra
2754EXAMPLE: example file2lplace; shows examples
2755"
2756{
2757  // format: from the usual string to letterplace
2759  // assume: file is a comma-sep list of polys
2760  // the vars are declared before
2761  // the file ends with ";"
2762  string t; int i;
2763  ideal I;
2764  list tst;
2765  while (s!="")
2766  {
2767    i = find(s,",");
2768    "i"; i;
2769    if (i==0)
2770    {
2771      i = find(s,";");
2772      if (i==0)
2773      {
2774        // no ; ??
2775         "no colon or semicolon found anymore";
2776         return(I);
2777      }
2778      // no "," but ";" on the i-th place
2779      t = s[1..i-1];
2780      s = "";
2781      "processing: "; t;
2782      tst = sent2lplace(t);
2783      tst;
2784      I = I, strList2poly(tst);
2785      return(I);
2786    }
2787    // here i !=0
2788    t = s[1..i-1];
2789    s = s[i+1..size(s)];
2790    "processing: "; t;
2791    tst = sent2lplace(t);
2792    tst;
2793    I = I, strList2poly(tst);
2794  }
2795  return(I);
2796}
2797example
2798{
2799  "EXAMPLE:"; echo = 2;
2800  ring r =0,(x,y,z,t),dp;
2801  def A = makeLetterplaceRing(4);
2802  setring A;
2803  string fn = "myfile";
2804  string s1 = "z*y*y*y - 3*y*z*x*y  + 3*y*y*z*y - y*x*y*z,";
2805  string s2 = "-2*y*x*y*z + y*y*z*z - z*z*y*y + 2*z*y*z*y,";
2806  string s3 = "z*y*x*t - 2*y*z*y*t + y*y*z*t - t*z*y*y + 2*t*y*z*y - t*x*y*z;";
2807  write(":w "+fn,s1);  write(":a "+fn,s2);   write(":a "+fn,s3);
2809  ideal I = file2lplace(fn);
2810  I;
2811}
2812
2813/* EXAMPLES AGAIN:
2814//static proc get_ls3nilp()
2815{
2816//first app of file2lplace
2817  ring r =0,(x,y,z,t),dp;
2818  int d = 10;
2819  def A = makeLetterplaceRing(d);
2820  setring A;
2821  ideal I = file2lplace("./ls3nilp.bg");
2822  // and now test the correctness: go back from lplace to strings
2823  lp2lstr(I,r);
2824  setring r;
2825  lst2str(@LN,1); // agree!
2826}
2827
2828//static proc doc_example()
2829{
2830  LIB "freegb.lib";
2831  ring r = 0,(x,y,z),dp;
2832  int d =4; // degree bound
2833  def R = makeLetterplaceRing(d);
2834  setring R;
2835  ideal I = x(1)*y(2) + y(1)*z(2), x(1)*x(2) + x(1)*y(2) - y(1)*x(2) - y(1)*y(2);
2836  option(redSB);option(redTail);
2837  ideal J = system("freegb",I,d,nvars(r));
2838  J;
2839  // visualization:
2840  lp2lstr(J,r); // export an object called @LN to the ring r
2841  setring r;  // change to the ring r
2842  lst2str(@LN,1); // output the strings
2843}
2844
2845*/
2846
2847// TODO:
2848// multiply two letterplace polynomials, lpMult: done
2849// reduction/ Normalform? needs kernel stuff
2850
2851proc lpMult(poly f, poly g)
2852"USAGE:  lpMult(f,g); f,g letterplace polynomials
2853RETURN:  poly
2854ASSUME: basering has a letterplace ring structure
2855PURPOSE: compute the letterplace form of f*g
2856EXAMPLE: example lpMult; shows examples
2857"
2858{
2859  if (lpAssumeViolation())
2860  {
2861    ERROR("Incomplete Letterplace structure on the basering!");
2862  }
2863  int sf = deg(f);
2864  int sg = deg(g);
2865  int uptodeg = attrib(basering, "uptodeg");
2866  if (sf+sg > uptodeg)
2867  {
2868    ERROR("degree bound violated by the product!");
2869  }
2870  //  if (sf>1) { sf = sf -1; }
2871  poly v = f*shiftPoly(g,sf);
2872  return(v);
2873}
2874example
2875{
2876  "EXAMPLE:"; echo = 2;
2877  // define a ring in letterplace form as follows:
2878  ring r = 0,(x(1),y(1),x(2),y(2),x(3),y(3),x(4),y(4)),dp;
2879  def R = setLetterplaceAttributes(r,4,2); // supply R with letterplace structure
2880  setring R;
2881  poly a = x(1)*y(2); poly b = y(1);
2882  lpMult(b,a);
2883  lpMult(a,b);
2884}
2885
2886proc lpPower(poly f, int n)
2887"USAGE:  lpPower(f,n); f letterplace polynomial, int n
2888RETURN:  poly
2889ASSUME: basering has a letterplace ring structure
2890PURPOSE: compute the letterplace form of f^n
2891EXAMPLE: example lpPower; shows examples
2892"
2893{
2894  if (n<0) { ERROR("the power must be a natural number!"); }
2895  if (n==0) { return(poly(1)); }
2896  if (n==1) { return(f); }
2897  int i;
2898  poly p = 1;
2899  for(i=1; i<= n; i++)
2900  {
2901    p = lpMult(p,f);
2902  }
2903  return(p);
2904}
2905example
2906{
2907  "EXAMPLE:"; echo = 2;
2908  // define a ring in letterplace form as follows:
2909  ring r = 0,(x(1),y(1),x(2),y(2),x(3),y(3),x(4),y(4)),dp;
2910  def R = setLetterplaceAttributes(r,4,2); // supply R with letterplace structure
2911  setring R;
2912  poly a = x(1)*y(2); poly b = y(1);
2913  lpPower(a,2);
2914  lpPower(b,4);
2915}
2916
2917// under development for Roberto
2918proc extractLinearPart(module M)
2919{
2920  /* returns vectors from a module whose max leadexp is 1 */
2921  /* does not take place nonlinearity into account yet */
2922  /* use rather kernel function isinV to get really nonlinear things */
2923  int i; int s = ncols(M);
2925  vector v; module Ret;
2926  for(i=1; i<=s; i++)
2927  {
2928    if ( isLinearVector(M[i]) )
2929    {
2930      Ret = Ret, M[i];
2931    }
2932  }
2933  Ret = simplify(Ret,2);
2934  return(Ret);
2935}
2936
2937// under development for Roberto
2938proc isLinearVector(vector v)
2939{
2940  /* returns true iff max leadexp is 1 */
2941  int i,j,k;
2942  intvec w;
2943  int s = size(v);
2944  poly p;
2946  for(i=1; i<=s; i++)
2947  {
2948    p = v[i];
2949    while (p != 0)
2950    {
2952      j = Max(w);
2953      if (j >=2)
2954      {
2957      }
2959    }
2960  }
2962}
2963
2964
2965// // the following is to determine a shift of a mono/poly from the
2966// // interface
2967
2968// proc whichshift(poly p, int numvars)
2969// {
2970// // numvars = number of vars of the orig free algebra
2971// // assume: we are in the letterplace ring
2972// // takes  monomial on the input
2975// if (v==0) { return(int(0)); }
2976// int sv = size(v);
2977// int i=1;
2978// while ( (v[i]==0) && (i<sv) ) { i++; }
2979// i = sv div i;
2980// return(i);
2981// }
2982
2983
2984// LIB "qhmoduli.lib";
2985// proc polyshift(poly p,  int numvars)
2986// {
2987//   poly q = p; int i = 0;
2988//   while (q!=0)
2989//   {
2990//     i = Max(i, whichshift(q,numvars));
2991//     q = q - lead(q);
2992//   }
2993//   return(q);
2994// }
2995
2996static proc lpAssumeViolation()
2997{
2998  // checks whether the global vars
2999  // uptodeg and lV are defined
3000  // returns Boolean : yes/no [for assume violation]
3001  def lpring = attrib(basering,"isLetterplaceRing");
3002  if ( typeof(lpring)!="int" )
3003  {
3004    //  if ( typeof(lpring)=="string" ) ??
3005    // basering is NOT lp Ring
3006
3007    return(1);
3008  }
3009  def uptodeg = attrib(basering,"uptodeg");
3010  if ( typeof(uptodeg)!="int" )
3011  {
3012    return(1);
3013  }
3014  def lV = attrib(basering,"lV");
3015  if ( typeof(lV)!="int" )
3016  {
3017    return(1);
3018  }
3019  //  int i = ( defined(uptodeg) && (defined(lV)) );
3020  //  return ( !i );
3021  return(0);
3022}
Note: See TracBrowser for help on using the repository browser.