# source:git/Singular/LIB/freegb.lib@a8a6b8

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