source: git/Singular/LIB/freegb.lib @ 86d5ad

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