source: git/Singular/LIB/freegb.lib @ 673fd0

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