source: git/Singular/LIB/freegb.lib @ 165e0e

spielwiese
Last change on this file since 165e0e was 165e0e, checked in by Karim Abou Zeid <karim23697@…>, 6 years ago
Update fpadim header, update fpadim and freegb test procedure
  • Property mode set to 100644
File size: 96.0 KB
Line 
1/////////////////////////////////////////////////////////////////////////////
2version="version freegb.lib 4.1.1.0 Dec_2017 "; // $Id$
3category="Noncommutative";
4info="
5LIBRARY: freegb.lib   Compute two-sided Groebner bases in free algebras via letterplace approach
6AUTHORS: Viktor Levandovskyy,     viktor.levandovskyy@math.rwth-aachen.de
7       Grischa Studzinski,      grischa.studzinski@math.rwth-aachen.de
8
9OVERVIEW: For the theory, see chapter 'Letterplace' in the @sc{Singular} Manual
10
11Support: Joint projects LE 2697/2-1 and KR 1907/3-1 of the Priority Programme SPP 1489:
12'Algorithmische und Experimentelle Methoden in Algebra, Geometrie und Zahlentheorie'
13of the German DFG
14and Project II.6 of the transregional collaborative research centre
15SFB-TRR 195 'Symbolic Tools in Mathematics and their Application' of the German DFG
16
17PROCEDURES:
18makeLetterplaceRing(d); creates a ring with d blocks of shifted original variables
19letplaceGBasis(I); computes two-sided Groebner basis of a letterplace ideal I up to a degree bound
20lpNF(f,I); two-sided normal form of f with respect to ideal I
21setLetterplaceAttributes(R,d,b); supplies ring R with the letterplace structure
22freeGBasis(L, n);  computes two-sided Groebner basis of an ideal, encoded via list L, up to degree n
23
24lpMult(f,g);    letterplace multiplication of letterplace polynomials
25shiftPoly(p,i); compute the i-th shift of letterplace polynomial p
26lpPower(f,n);   natural power of a letterplace polynomial
27lieBracket(a,b[, N]);  compute Lie bracket ab-ba of two letterplace polynomials
28
29lp2lstr(K, s);  convert a letterplace ideal into a list of modules
30lst2str(L[, n]); convert a list (of modules) into polynomials in free algebra via strings
31mod2str(M[, n]); convert a module into a polynomial in free algebra via strings
32vct2str(M[, n]);   convert a vector into a word in free algebra
33
34serreRelations(A,z); compute the homogeneous part of Serre's relations associated to a generalized Cartan matrix A
35fullSerreRelations(A,N,C,P,d); compute the ideal of all Serre's relations associated to a generalized Cartan matrix A
36isVar(p);              check whether p is a power of a single variable
37ademRelations(i,j);    compute the ideal of Adem relations for i<2j in char 0
38
39SEE ALSO: fpadim_lib, LETTERPLACE
40";
41
42// this library computes two-sided GB of an ideal
43// in a free associative algebra
44
45// a monomial is encoded via a vector V
46// where V[1] = coefficient
47// V[1+i] = the corresponding symbol
48
49LIB "qhmoduli.lib"; // for Max
50LIB "bfun.lib"; // for inForm
51
52proc tstfreegb()
53{
54    /* tests all procs for consistency */
55  /* adding the new proc, add it here */
56  example makeLetterplaceRing;
57  example letplaceGBasis;
58  example 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 = 2; // temporary until makeLetterplaceRing4() is fixed
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  if (deg(a) + i > uptodeg)
2524  {
2525    ERROR("degree bound violated by the shift!");
2526  }
2527  return(system("stest",a,i,uptodeg,lV));
2528}
2529example
2530{
2531  "EXAMPLE:"; echo = 2;
2532  ring r = 0,(x,y,z),dp;
2533  int uptodeg = 5; int lV = 3;
2534  def R = makeLetterplaceRing(uptodeg);
2535  setring R;
2536  poly f = x(1)*z(2)*y(3) - 2*z(1)*y(2) + 3*x(1);
2537  shiftPoly(f,1);
2538  shiftPoly(f,2);
2539}
2540
2541
2542static proc mmLiebr(poly a, poly b)
2543{
2544  // a,b, monomials
2545  a = lead(a);
2546  b = lead(b);
2547  int sa = deg(a);
2548  int sb = deg(b);
2549  poly v = a*shiftPoly(b,sa) - b*shiftPoly(a,sb);
2550  return(v);
2551}
2552
2553static proc test_shift()
2554{
2555  LIB "freegb.lib";
2556  ring r = 0,(a,b),dp;
2557  int d =5;
2558  def R = makeLetterplaceRing(d);
2559  setring R;
2560  int uptodeg = d;
2561  int lV = 2;
2562  def R = setLetterplaceAttributes(r,uptodeg,2); // supply R with letterplace structure
2563  setring R;
2564  poly p = mmLiebr(a(1),b(1));
2565  poly p = lieBracket(a(1),b(1));
2566}
2567
2568proc serreRelations(intmat A, int zu)
2569"USAGE:  serreRelations(A,z); A an intmat, z an int
2570RETURN:  ideal
2571ASSUME: basering has a letterplace ring structure and
2572@*          A is a generalized Cartan matrix with integer entries
2573PURPOSE: compute the ideal of Serre's relations associated to A
2574EXAMPLE: example serreRelations; shows examples
2575"
2576{
2577  // zu = 1 -> with commutators [f_i,f_j]; zu == 0 without them
2578  // suppose that A is cartan matrix
2579  // then Serre's relations are
2580  // (ad f_j)^{1-A_{ij}} ( f_i)
2581  int ppl = printlevel-voice+2;
2582  int n = ncols(A); // hence n variables
2583  int i,j,k,el;
2584  poly p,q;
2585  ideal I;
2586  for (i=1; i<=n; i++)
2587  {
2588    for (j=1; j<=n; j++)
2589    {
2590      el = 1 - A[i,j];
2591      //     printf("i:%s, j: %s, l: %s",i,j,l);
2592      dbprint(ppl,"i, j, l: ",i,j,el);
2593      //      if ((i!=j) && (l >0))
2594      //      if ( (i!=j) &&  ( ((zu ==0) &&  (l >=2)) || ((zu ==1) &&  (l >=1)) ) )
2595      if ((i!=j) && (el >0))
2596      {
2597        q = lieBracket(var(j),var(i));
2598        dbprint(ppl,"first bracket: ",q);
2599        //        if (l >=2)
2600        //        {
2601          for (k=1; k<=el-1; k++)
2602          {
2603            q = lieBracket(var(j),q);
2604            dbprint(ppl,"further bracket:",q);
2605          }
2606          //        }
2607      }
2608      if (q!=0) { I = I,q; q=0;}
2609    }
2610  }
2611  I = simplify(I,2);
2612  return(I);
2613}
2614example
2615{
2616  "EXAMPLE:"; echo = 2;
2617  intmat A[3][3] =
2618    2, -1, 0,
2619    -1, 2, -3,
2620    0, -1, 2; // G^1_2 Cartan matrix
2621  ring r = 0,(f1,f2,f3),dp;
2622  int uptodeg = 5;
2623  def R = makeLetterplaceRing(uptodeg);
2624  setring R;
2625  ideal I = serreRelations(A,1); I = simplify(I,1+2+8);
2626  I;
2627}
2628
2629/* setup for older example:
2630  intmat A[2][2] = 2, -1, -1, 2; // sl_3 == A_2
2631  ring r = 0,(f1,f2),dp;
2632  int uptodeg = 5; int lV = 2;
2633*/
2634
2635proc fullSerreRelations(intmat A, ideal rNegative, ideal rCartan, ideal rPositive, int uptodeg)
2636"USAGE:  fullSerreRelations(A,N,C,P,d); A an intmat, N,C,P ideals, d an int
2637RETURN:  ring (and ideal)
2638PURPOSE: compute the inhomogeneous Serre's relations associated to A in given
2639@*       variable names
2640ASSUME: three ideals in the input are of the same sizes and contain merely
2641@* variables which are interpreted as follows: N resp. P stand for negative
2642@* resp. positive roots, C stand for Cartan elements. d is the degree bound for
2643@* letterplace ring, which will be returned.
2644@* The matrix A is a generalized Cartan matrix with integer entries
2645@* The result is the ideal called 'fsRel' in the returned ring.
2646EXAMPLE: example fullSerreRelations; shows examples
2647"
2648{
2649  /* SerreRels on rNeg and rPos plus Cartans etc. */
2650  int ppl = printlevel -voice+2;
2651  /* ideals must be written in variables: assume each term is of degree 1 */
2652  int i,j,k;
2653  int N = nvars(basering);
2654  def save = basering;
2655  int comFlag = 0;
2656  /* assume:  (size(rNegative) == size(rPositive)) */
2657  /* assume:  (size(rNegative) == size(rCartan)) i.e. nonsimple Cartans */
2658  if ( (size(rNegative) != size(rPositive)) || (size(rNegative) != size(rCartan)) )
2659  {
2660    ERROR("All input ideals must be of the same size");
2661  }
2662
2663//   if (size(rNegative) != size(rPositive))
2664//   {
2665//     ERROR("The 1st and the 3rd input ideals must be of the same size");
2666//   }
2667
2668  /* assume:  2*size(rNegative) + size(rCartan) >= nvars(basering) */
2669  i = 2*size(rNegative) + size(rCartan);
2670  if (i>N)
2671  {
2672    string s1="The total number of elements in input ideals";
2673    string s2="must not exceed the dimension of the ground ring";
2674    ERROR(s1+s2);
2675  }
2676  if (i < N)
2677  {
2678    comFlag = N-i; // so many elements will commute
2679    "Warning: some elements will be treated as mutually commuting";
2680  }
2681  /* extract varnames from input ideals */
2682  intvec iNeg = varIdeal2intvec(rNegative);
2683  intvec iCartan = varIdeal2intvec(rCartan);
2684  intvec iPos = varIdeal2intvec(rPositive);
2685  /* for each vector in rNeg and rPositive, go into the corr. ring and create SerreRels */
2686  /* rNegative: */
2687  list L = ringlist(save);
2688  def LPsave = makeLetterplaceRing2(uptodeg); setring save;
2689  list LNEG = L; list tmp;
2690  /* L[1] field as is; L[2] vars: a subset; L[3] ordering: dp, L[4] as is */
2691  for (i=1; i<=size(iNeg); i++)
2692  {
2693    tmp[i] = string(var(iNeg[i]));
2694  }
2695  LNEG[2] = tmp; LNEG[3] = list(list("dp",intvec(1:size(iNeg))), list("C",0));
2696  def RNEG = ring(LNEG); setring RNEG;
2697  def RRNEG = makeLetterplaceRing2(uptodeg);
2698  setring RRNEG;
2699  ideal I = serreRelations(A,1); I = simplify(I,1+2+8);
2700  setring LPsave;
2701  ideal srNeg = imap(RRNEG,I);
2702  dbprint(ppl,"0-1 ideal of negative relations is ready");
2703  dbprint(ppl-1,srNeg);
2704  setring save; kill L,tmp,RRNEG,RNEG, LNEG;
2705  /* rPositive: */
2706  list L = ringlist(save);
2707  list LPOS = L; list tmp;
2708  /* L[1] field as is; L[2] vars: a subset; L[3] ordering: dp, L[4] as is */
2709  for (i=1; i<=size(iPos); i++)
2710  {
2711    tmp[i] = string(var(iPos[i]));
2712  }
2713  LPOS[2] = tmp; LPOS[3] = list(list("dp",intvec(1:size(iPos))), list("C",0));
2714  def RPOS = ring(LPOS); setring RPOS;
2715  def RRPOS = makeLetterplaceRing2(uptodeg);
2716  setring RRPOS;
2717  ideal I = serreRelations(A,1); I = simplify(I,1+2+8);
2718  setring LPsave;
2719  ideal srPos = imap(RRPOS,I);
2720  dbprint(ppl,"0-2 ideal of positive relations is ready");
2721  dbprint(ppl-1,srPos);
2722  setring save; kill L,tmp,RRPOS,RPOS, LPOS;
2723  string sMap = "ideal Mmap =";
2724  for (i=1; i<=nvars(save); i++)
2725  {
2726    sMap = sMap + string(var(i)) +"(1),";
2727  }
2728  sMap[size(sMap)] = ";";
2729  /* cartans: h_j h_i = h_i h_j */
2730  setring LPsave;
2731  ideal ComCartan;
2732  for (i=1; i<size(iCartan); i++)
2733  {
2734    for (j=i+1; j<=size(iCartan); j++)
2735    {
2736      ComCartan =  ComCartan + lieBracket(var(iCartan[j]),var(iCartan[i]));
2737    }
2738  }
2739  ComCartan = simplify(ComCartan,1+2+8);
2740  execute(sMap); // defines an ideal Mmap
2741  map F = save, Mmap;
2742  dbprint(ppl,"1. commuting Cartans: ");
2743  dbprint(ppl-1,ComCartan);
2744  /* [e_i, f_j] =0 if i<>j */
2745  ideal ComPosNeg; // assume: #Neg=#Pos
2746  for (i=1; i<size(iPos); i++)
2747  {
2748    for (j=1; j<=size(iPos); j++)
2749    {
2750      if (j !=i)
2751      {
2752        ComPosNeg =  ComPosNeg + lieBracket(var(iPos[i]),var(iNeg[j]));
2753        ComPosNeg =  ComPosNeg + lieBracket(var(iPos[j]),var(iNeg[i]));
2754      }
2755    }
2756  }
2757  ComPosNeg = simplify(ComPosNeg,1+2+8);
2758  dbprint(ppl,"2. commuting Positive and Negative:");
2759  dbprint(ppl-1,ComPosNeg);
2760  /* [e_i, f_i] = h_i */
2761  poly tempo;
2762  for (i=1; i<=size(iCartan); i++)
2763  {
2764    tempo = lieBracket(var(iPos[i]),var(iNeg[i])) - var(iCartan[i]);
2765    ComPosNeg =  ComPosNeg + tempo;
2766  }
2767  //  ComPosNeg = simplify(ComPosNeg,1+2+8);
2768  dbprint(ppl,"3. added sl2 triples [e_i,f_i]=h_i");
2769  dbprint(ppl-1,ComPosNeg);
2770
2771  /* [h_i, e_j] = A_ij e_j */
2772  /* [h_i, f_j] = -A_ij f_j */
2773  ideal ActCartan; // assume: #Neg=#Pos
2774  for (i=1; i<=size(iCartan); i++)
2775  {
2776    for (j=1; j<=size(iCartan); j++)
2777    {
2778      tempo = lieBracket(var(iCartan[i]),var(iPos[j])) - A[i,j]*var(iPos[j]);
2779      ActCartan = ActCartan + tempo;
2780      tempo = lieBracket(var(iCartan[i]),var(iNeg[j])) + A[i,j]*var(iNeg[j]);
2781      ActCartan = ActCartan + tempo;
2782    }
2783  }
2784  ActCartan = simplify(ActCartan,1+2+8);
2785  dbprint(ppl,"4. actions of Cartan:");
2786  dbprint(ppl-1, ActCartan);
2787
2788  /* final part: prepare the output */
2789  setring LPsave;
2790  ideal fsRel = srNeg, srPos, ComPosNeg, ComCartan, ActCartan;
2791  export fsRel;
2792  setring save;
2793  return(LPsave);
2794}
2795example
2796{
2797  "EXAMPLE:"; echo = 2;
2798  intmat A[2][2] =
2799    2, -1,
2800    -1, 2; // A_2 = sl_3 Cartan matrix
2801  ring r = 0,(f1,f2,h1,h2,e1,e2),dp;
2802  ideal negroots = f1,f2; ideal cartans = h1,h2; ideal posroots = e1,e2;
2803  int uptodeg = 5;
2804  def RS = fullSerreRelations(A,negroots,cartans,posroots,uptodeg);
2805  setring RS; fsRel;
2806}
2807
2808static proc varIdeal2intvec(ideal I)
2809{
2810  // used in SerreRelations
2811  /* assume1:  input ideal is a list of variables of the ground ring */
2812  int i,j; intvec V;
2813  for (i=1; i<= size(I); i++)
2814  {
2815    j = univariate(I[i]);
2816    if (j<=0)
2817    {
2818      ERROR("input ideal must contain only variables");
2819    }
2820    V[i] = j;
2821  }
2822  dbprint(printlevel-voice+2,V);
2823  /* now we make a smaller list of non-repeating entries */
2824  ideal iW = simplify(ideal(V),2+4); // no zeros, no repetitions
2825  if (size(iW) < size(V))
2826  {
2827    /* extract intvec from iW */
2828    intvec inW;
2829    for(j=1; j<=size(iW); j++)
2830    {
2831      inW[j] = int(leadcoef(iW[j]));
2832    }
2833    return(inW);
2834  }
2835  return(V);
2836}
2837example
2838{
2839  "EXAMPLE:"; echo = 2;
2840  ring r = 0,(x,y,z),dp;
2841  ideal I = x,z;
2842  varIdeal2intvec(I);
2843  varIdeal2intvec(ideal(x2,y^3,x+1));
2844  varIdeal2intvec(ideal(x*y,y,x+1));
2845}
2846
2847proc lp2lstr(ideal K, def save)
2848"USAGE:  lp2lstr(K,s); K an ideal, s a ring name
2849RETURN:  nothing (exports object @LN into the ring named s)
2850ASSUME: basering has a letterplace ring structure
2851PURPOSE: converts letterplace ideal to list of modules
2852NOTE: useful as preprocessing to 'lst2str'
2853EXAMPLE: example lp2lstr; shows examples
2854"
2855{
2856  def @R = basering;
2857  string err;
2858  int s = nvars(save);
2859  int i,j,k;
2860    // K contains vars x(1),...z(1) = images of originals
2861  // 5. go back to orig vars, produce strings/modules
2862  int sk = size(K);
2863  int sp, sx, a, b;
2864  intvec x;
2865  poly p,q;
2866  poly pn;
2867  // vars in 'save'
2868  setring save;
2869  module N;
2870  list LN;
2871  vector V;
2872  poly pn;
2873  // test and skip exponents >=2
2874  setring @R;
2875  for(i=1; i<=sk; i++)
2876  {
2877    p  = K[i];
2878    while (p!=0)
2879    {
2880      q  = lead(p);
2881      //      "processing q:";q;
2882      x  = leadexp(q);
2883      sx = size(x);
2884      for(k=1; k<=sx; k++)
2885      {
2886        if ( x[k] >= 2 )
2887        {
2888          err = "skip: the value x[k] is " + string(x[k]);
2889          dbprint(ppl,err);
2890          //            return(0);
2891          K[i] = 0;
2892          p    = 0;
2893          q    = 0;
2894          break;
2895        }
2896      }
2897      p  = p - q;
2898    }
2899  }
2900  K  = simplify(K,2);
2901  sk = size(K);
2902  for(i=1; i<=sk; i++)
2903  {
2904    //    setring save;
2905    //    V  = 0;
2906    setring @R;
2907    p  = K[i];
2908    while (p!=0)
2909    {
2910      q  = lead(p);
2911      err =  "processing q:" + string(q);
2912      dbprint(ppl,err);
2913      x  = leadexp(q);
2914      sx = size(x);
2915      pn = leadcoef(q);
2916      setring save;
2917      pn = imap(@R,pn);
2918      V  = V + leadcoef(pn)*gen(1);
2919      for(k=1; k<=sx; k++)
2920      {
2921        if (x[k] ==1)
2922        {
2923          a = k div s; // block number=a+1, a!=0
2924          b = k % s; // remainder
2925          //          printf("a: %s, b: %s",a,b);
2926          if (b == 0)
2927          {
2928            // that is it's the last var in the block
2929            b = s;
2930            a = a-1;
2931          }
2932          V = V + var(b)*gen(a+2);
2933        }
2934      }
2935      err = "V: " + string(V);
2936      dbprint(ppl,err);
2937      //      printf("V: %s", string(V));
2938      N = N,V;
2939      V  = 0;
2940      setring @R;
2941      p  = p - q;
2942      pn = 0;
2943    }
2944    setring save;
2945    LN[i] = simplify(N,2);
2946    N     = 0;
2947  }
2948  setring save;
2949  list @LN = LN;
2950  export @LN;
2951  //  return(LN);
2952}
2953example
2954{
2955  "EXAMPLE:"; echo = 2;
2956  intmat A[2][2] = 2, -1, -1, 2; // sl_3 == A_2
2957  ring r = 0,(f1,f2),dp;
2958  def R = makeLetterplaceRing(3);
2959  setring R;
2960  ideal I = serreRelations(A,1);
2961  lp2lstr(I,r);
2962  setring r;
2963  lst2str(@LN,1);
2964}
2965
2966static proc strList2poly(list L)
2967{
2968  //  list L comes from sent2lplace (which takes a polynomial as input)
2969  // each entry of L is a sublist with the coef on the last place
2970  int s = size(L); int t;
2971  int i,j;
2972  list M;
2973  poly p,q;
2974  string Q;
2975  for(i=1; i<=s; i++)
2976  {
2977    M = L[i];
2978    t = size(M);
2979    //    q = M[t]; // a constant
2980    Q = string(M[t]);
2981    for(j=1; j<t; j++)
2982    {
2983      //      q = q*M[j];
2984      Q = Q+"*"+string(M[j]);
2985    }
2986    execute("q="+Q+";");
2987    //    q;
2988    p = p + q;
2989  }
2990  kill Q;
2991  return(p);
2992}
2993example
2994{
2995  "EXAMPLE:"; echo = 2;
2996  ring r =0,(x,y,z,t),Dp;
2997  def A = makeLetterplaceRing(4);
2998  setring A;
2999  string t = "-2*y*z*y*z + y*t*z*z - z*x*x*y  + 2*z*y*z*y";
3000  list L = sent2lplace(t);
3001  L;
3002  poly p = strList2poly(L);
3003  p;
3004}
3005
3006static proc file2lplace(string fname)
3007"USAGE:  file2lplace(fnm);  fnm a string
3008RETURN:  ideal
3009PURPOSE: convert the contents of the file fnm into ideal of polynomials in free algebra
3010EXAMPLE: example file2lplace; shows examples
3011"
3012{
3013  // format: from the usual string to letterplace
3014  string s = read(fname);
3015  // assume: file is a comma-sep list of polys
3016  // the vars are declared before
3017  // the file ends with ";"
3018  string t; int i;
3019  ideal I;
3020  list tst;
3021  while (s!="")
3022  {
3023    i = find(s,",");
3024    "i"; i;
3025    if (i==0)
3026    {
3027      i = find(s,";");
3028      if (i==0)
3029      {
3030        // no ; ??
3031         "no colon or semicolon found anymore";
3032         return(I);
3033      }
3034      // no "," but ";" on the i-th place
3035      t = s[1..i-1];
3036      s = "";
3037      "processing: "; t;
3038      tst = sent2lplace(t);
3039      tst;
3040      I = I, strList2poly(tst);
3041      return(I);
3042    }
3043    // here i !=0
3044    t = s[1..i-1];
3045    s = s[i+1..size(s)];
3046    "processing: "; t;
3047    tst = sent2lplace(t);
3048    tst;
3049    I = I, strList2poly(tst);
3050  }
3051  return(I);
3052}
3053example
3054{
3055  "EXAMPLE:"; echo = 2;
3056  ring r =0,(x,y,z,t),dp;
3057  def A = makeLetterplaceRing(4);
3058  setring A;
3059  string fn = "myfile";
3060  string s1 = "z*y*y*y - 3*y*z*x*y  + 3*y*y*z*y - y*x*y*z,";
3061  string s2 = "-2*y*x*y*z + y*y*z*z - z*z*y*y + 2*z*y*z*y,";
3062  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;";
3063  write(":w "+fn,s1);  write(":a "+fn,s2);   write(":a "+fn,s3);
3064  read(fn);
3065  ideal I = file2lplace(fn);
3066  I;
3067}
3068
3069/* EXAMPLES AGAIN:
3070//static proc get_ls3nilp()
3071{
3072//first app of file2lplace
3073  ring r =0,(x,y,z,t),dp;
3074  int d = 10;
3075  def A = makeLetterplaceRing(d);
3076  setring A;
3077  ideal I = file2lplace("./ls3nilp.bg");
3078  // and now test the correctness: go back from lplace to strings
3079  lp2lstr(I,r);
3080  setring r;
3081  lst2str(@LN,1); // agree!
3082}
3083
3084//static proc doc_example()
3085{
3086  LIB "freegb.lib";
3087  ring r = 0,(x,y,z),dp;
3088  int d =4; // degree bound
3089  def R = makeLetterplaceRing(d);
3090  setring R;
3091  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);
3092  option(redSB);option(redTail);
3093  ideal J = system("freegb",I,d,nvars(r));
3094  J;
3095  // visualization:
3096  lp2lstr(J,r); // export an object called @LN to the ring r
3097  setring r;  // change to the ring r
3098  lst2str(@LN,1); // output the strings
3099}
3100
3101*/
3102
3103static proc lpMultX(poly f, poly g)
3104{
3105  /* multiplies two polys in a very general setting correctly */
3106  /* alternative to lpMult, possibly better at non-positive orderings */
3107
3108  if (lpAssumeViolation())
3109  {
3110    ERROR("Incomplete Letterplace structure on the basering!");
3111  }
3112  // decompose f,g into graded pieces with inForm: need dmodapp.lib
3113  int b = attrib(basering,"lV");  // the length of the block
3114  intvec w; // inherit the graded on the oridinal ring
3115  int i;
3116  for(i=1; i<=b; i++)
3117  {
3118    w[i] = deg(var(i));
3119  }
3120  intvec v = w;
3121  for(i=1; i< attrib(basering,"uptodeg"); i++)
3122  {
3123    v = v,w;
3124  }
3125  w = v;
3126  poly p,q,s, result;
3127  s = g;
3128  while (f!=0)
3129  {
3130    p = inForm(f,w)[1];
3131    f = f - p;
3132    s = g;
3133    while (s!=0)
3134    {
3135      q = inForm(s,w)[1];
3136      s = s - q;
3137      result = result + lpMult(p,q);
3138    }
3139  }
3140  // shrinking
3141  //  result;
3142  return( system("shrinktest",result,attrib(basering, "lV")) );
3143}
3144example
3145{
3146  "EXAMPLE:"; echo = 2;
3147  // define a ring in letterplace form as follows:
3148  ring r = 0,(x(1),y(1),x(2),y(2),x(3),y(3),x(4),y(4)),dp;
3149  def R = setLetterplaceAttributes(r,4,2); // supply R with letterplace structure
3150  setring R;
3151  poly a = x(1)*y(2)+x(1)+y(1); poly b = y(1)+3;
3152  lpMultX(b,a);
3153  lpMultX(a,b);
3154}
3155
3156// multiply two letterplace polynomials, lpMult: done
3157// reduction/ Normalform? needs kernel stuff
3158
3159
3160proc lpMult(poly f, poly g)
3161"USAGE:  lpMult(f,g); f,g letterplace polynomials
3162RETURN:  poly
3163ASSUME: basering has a letterplace ring structure
3164PURPOSE: compute the letterplace form of f*g
3165EXAMPLE: example lpMult; shows examples
3166"
3167{
3168
3169  // changelog:
3170  // VL oct 2010: deg -> deg(_,w) for the length
3171  // shrink the result => don't need to decompose polys
3172  // since the shift is big enough
3173
3174  // indeed it's better to have that
3175  // ASSUME: both f and g are quasi-homogeneous
3176
3177  if (lpAssumeViolation())
3178  {
3179    ERROR("Incomplete Letterplace structure on the basering!");
3180  }
3181  intvec w = 1:nvars(basering);
3182  int sf = deg(f,w); // VL Oct 2010: we need rather length than degree
3183  int sg = deg(g,w); // esp. in the case of weighted ordering
3184  int uptodeg = attrib(basering, "uptodeg");
3185  if (sf+sg > uptodeg)
3186  {
3187    ERROR("degree bound violated by the product!");
3188  }
3189  //  if (sf>1) { sf = sf -1; }
3190  poly v = f*shiftPoly(g,sf);
3191  // bug, reported by Simon King: in nonhomog case [solved]
3192  // we need to shrink
3193  return( system("shrinktest",v,attrib(basering, "lV")) );
3194}
3195example
3196{
3197  "EXAMPLE:"; echo = 2;
3198  // define a ring in letterplace form as follows:
3199  ring r = 0,(x(1),y(1),x(2),y(2),x(3),y(3),x(4),y(4)),dp;
3200  def R = setLetterplaceAttributes(r,4,2); // supply R with letterplace structure
3201  setring R;
3202  poly a = x(1)*y(2)+x(1)+y(1); poly b = y(1)+3;
3203  lpMult(b,a);
3204  lpMult(a,b);
3205}
3206
3207proc lpPower(poly f, int n)
3208"USAGE:  lpPower(f,n); f letterplace polynomial, int n
3209RETURN:  poly
3210ASSUME: basering has a letterplace ring structure
3211PURPOSE: compute the letterplace form of f^n
3212EXAMPLE: example lpPower; shows examples
3213"
3214{
3215  if (n<0) { ERROR("the power must be a natural number!"); }
3216  if (n==0) { return(poly(1)); }
3217  if (n==1) { return(f); }
3218  int i;
3219  poly p = 1;
3220  for(i=1; i<= n; i++)
3221  {
3222    p = lpMult(p,f);
3223  }
3224  return(p);
3225}
3226example
3227{
3228  "EXAMPLE:"; echo = 2;
3229  // define a ring in letterplace form as follows:
3230  ring r = 0,(x(1),y(1),x(2),y(2),x(3),y(3),x(4),y(4)),dp;
3231  def R = setLetterplaceAttributes(r,4,2); // supply R with letterplace structure
3232  setring R;
3233  poly a = x(1)*y(2) + y(1); poly b = y(1) - 1;
3234  lpPower(a,2);
3235  lpPower(b,4);
3236}
3237
3238// new: lp normal from by using shift-invariant data by Grischa Studzinski
3239
3240/////////////////////////////////////////////////////////
3241//   ASSUMPTIONS: every polynomial is an element of V',
3242//@* else there wouldn't be an dvec representation
3243
3244//Main procedure for the user
3245
3246proc lpNF(poly p, ideal G)
3247"USAGE: lpNF(p,G); f letterplace polynomial, ideal I
3248RETURN: poly
3249PURPOSE: computation of the normal form of p with respect to G
3250ASSUME: p is a Letterplace polynomial, G is a set Letterplace polynomials,
3251being a Letterplace Groebner basis (no check for this will be done)
3252NOTE: Strategy: take the smallest monomial wrt ordering for reduction
3253-     For homogenous ideals the shift does not matter
3254-     For non-homogenous ideals the first shift will be the smallest monomial
3255EXAMPLE: example lpNF; shows examples
3256"
3257{if ((p==0) || (size(G) == 0)){return(p);}
3258 checkAssumptions(p,G);
3259 G = sort(G)[1];
3260 list L = makeDVecI(G);
3261 return(normalize(lpNormalForm2(p,G,L)));
3262}
3263example
3264{
3265  "EXAMPLE:"; echo = 2;
3266ring r = 0,(x,y,z),dp;
3267int d =5; // degree
3268def R = makeLetterplaceRing(d);
3269setring R;
3270ideal 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);
3271ideal J = letplaceGBasis(I); // compute a Letterplace Groebner basis
3272poly 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);
3273poly q = z(1)*x(2)*z(3)*y(4)*z(5) - y(1)*z(2)*x(3)*y(4)*z(5);
3274lpNF(p,J);
3275lpNF(q,J);
3276}
3277
3278//procedures to convert monomials into the DVec representation, all static
3279////////////////////////////////////////////////////////
3280
3281
3282static proc getExpVecs(ideal G)
3283"USUAGE: getExpVecs(G);
3284RETURN: list of intvecs
3285PURPOSE: convert G into a list of intvecs, corresponding to the exponent vector
3286 of the leading monomials of G
3287"
3288{int i; list L;
3289 for (i = 1; i <= size(G); i++) {L[i] = leadexp(G[i]); }
3290 return(L);
3291}
3292
3293static proc delSupZero(intvec I)
3294"USUAGE:delSupZero(I);
3295RETURN: intvec
3296PURPOSE: Deletes superfluous zero blocks of an exponent vector
3297ASSUME: Intvec is an exponent vector of a letterplace monomial contained in V'
3298"
3299{if (I==intvec(0)) {return(intvec(0));}
3300 int j,k,l;
3301 int n = attrib(basering,"lV"); int d = attrib(basering,"uptodeg");
3302 intvec w; j = 1;
3303 while (j <= d)
3304 {w = I[1..n];
3305  if (w<>intvec(0)){break;}
3306   else {I = I[(n+1)..(n*d)]; d = d-1; j++;}
3307 }
3308 for (j = 1; j <= d; j++)
3309  {l=(j-1)*n+1; k= j*n;
3310   w = I[l..k];
3311   if (w==intvec(0)){w = I[1..(l-1)]; return(w);}//if a zero block is found there are only zero blocks left,
3312                                                 //otherwise there would be a hole in the monomial
3313                                                 // shrink should take care that this will not happen
3314  }
3315 return(I);
3316}
3317
3318static proc delSupZeroList(list L)
3319"USUAGE:delSupZeroList(L); L a list, containing intvecs
3320RETURN: list, containing intvecs
3321PURPOSE: Deletes all superfluous zero blocks for a list of exponent vectors
3322ASSUME: All intvecs are exponent vectors of letterplace monomials contained in V'
3323"
3324{int i;
3325 for (i = size(L); 0 < i; i--){L[i] = delSupZero(L[i]);}
3326 return(L);
3327}
3328
3329
3330static proc makeDVec(intvec V)
3331"USUAGE:makeDVec(V);
3332RETURN: intvec
3333PURPOSE: Converts an modified exponent vector into an Dvec
3334NOTE: Superfluos zero blocks must have been deleted befor using this procedure
3335"
3336{int i,j,k,r1,r2; intvec D;
3337 int n = attrib(basering,"lV");
3338 k = size(V) div n; r1 = 0; r2 = 0;
3339 for (i=1; i<= k; i++)
3340  {for (j=(1+((i-1)*n)); j <= (i*n); j++)
3341   {if (V[j]>0){r2 = j - ((i-1)*n); j = (j mod n); break;}
3342   }
3343   D[size(D)+1] = r1+r2;
3344   if (j == 0) {r1 = 0;} else{r1= n-j;}
3345  }
3346 D = D[2..size(D)];
3347 return(D);
3348}
3349
3350static proc makeDVecL(list L)
3351"USUAGE:makeDVecL(L); L, a list containing intvecs
3352RETURN: list, containing intvecs
3353ASSUME:
3354"
3355{int i; list R;
3356 for (i=1; i <= size(L); i++) {R[i] = makeDVec(L[i]);}
3357 return(R);
3358}
3359
3360static proc makeDVecI(ideal G)
3361"USUAGE:makeDVecI(G);
3362RETURN:list, containing intvecs
3363PURPOSE:computing the DVec representation for lead(G)
3364ASSUME:
3365"
3366{list L = delSupZeroList(getExpVecs(G));
3367 return(makeDVecL(L));
3368}
3369
3370
3371//procedures, which are dealing with the DVec representation, all static
3372
3373static proc dShiftDiv(intvec V, intvec W)
3374"USUAGE: dShiftDiv(V,W);
3375RETURN: a list,containing integers, or -1, if no shift of W divides V
3376PURPOSE: find all possible shifts s, such that s.W|V
3377ASSUME: V,W are DVecs of monomials contained in V'
3378"
3379{if(size(V)<size(W)){return(list(-1));}
3380
3381 int i,j,r; intvec T; list R;
3382 int n = attrib(basering,"lV");
3383 int k = size(V) - size(W) + 1;
3384 if (intvec(V[1..size(W)])-W == 0){R[1]=0;}
3385 for (i =2; i <=k; i++)
3386 {r = 0; kill T; intvec T;
3387  for (j =1; j <= i; j++) {r = r + V[j];}
3388  //if (i==1) {T[1] = r-(i-1)*n;} else
3389  T[1] = r-(i-1)*n; if (size(W)>1) {T[2..size(W)] = V[(i+1)..(size(W)+i-1)];}
3390  if (T-W == 0) {R[size(R)+1] = i-1;}
3391 }
3392 if (size(R)>0) {return(R);}
3393 else {return(list(-1));}
3394}
3395
3396//the first normal form procedure, if a user want not to presort the ideal, just make it not static
3397
3398static proc lpNormalForm1(poly p, ideal G, list L)
3399"USUAGE:lpNormalForm1(p,G);
3400RETURN:poly
3401PURPOSE:computation of the normalform of p w.r.t. G
3402ASSUME: p is a Letterplace polynomial, G is a set of Letterplace polynomials
3403NOTE: Taking the first possible reduction
3404"
3405{
3406 if (deg(p) <1) {return(p);}
3407 else
3408  {
3409  int i; int s;
3410  intvec V = makeDVec(delSupZero(leadexp(p)));
3411  for (i = 1; i <= size(L); i++)
3412  {s = dShiftDiv(V, L[i])[1];
3413   if (s <> -1)
3414   {p = lpReduce(p,G[i],s);
3415    p = lpNormalForm1(p,G,L);
3416    break;
3417   }
3418  }
3419  p = p[1] + lpNormalForm1(p-p[1],G,L);
3420  return(p);
3421 }
3422}
3423
3424
3425// new VL; called from lpNF
3426static proc lpNormalForm2(poly pp, ideal G, list L)
3427"USUAGE:lpNormalForm2(p,G);
3428RETURN:poly
3429PURPOSE:computation of the normal form of p w.r.t. G
3430ASSUME: p is a Letterplace polynomial, G is a set of Letterplace polynomials
3431NOTE: Taking the first possible reduction
3432"
3433{
3434 poly one = 1;
3435 if ( (pp == 0) || (leadmonom(pp) == one) ) { return(pp); }
3436 poly p = pp; poly q;
3437 int i; int s; intvec V;
3438 while ( (p != 0) && (leadmonom(p) != one) )
3439 {
3440   //"entered while with p="; p;
3441   V = makeDVec(delSupZero(leadexp(p)));
3442   i = 0;
3443   s = -1;
3444   //"look for divisor";
3445   while ( (s == -1) && (i<size(L)) )
3446   {
3447     i = i+1;
3448     s = dShiftDiv(V, L[i])[1];
3449   }
3450 // now, out of here: either i=size(L) and s==-1 => no reduction
3451 // otherwise: i<=size(L) and s!= -1 => reduction
3452    //"out of divisor search: s="; s; "i="; i;
3453    if (s != -1)
3454    {
3455    //"start reducing with G[i]:";
3456      p = lpReduce(p,G[i],s); // lm-reduction
3457      //"reduced to p="; p;
3458    }
3459    else
3460    {
3461      // ie no lm-reduction possible; proceed with the tail reduction
3462      q = p-lead(p);
3463      p = lead(p);
3464      if (q!=0)
3465      {
3466        p = p + lpNormalForm2(q,G,L);
3467      }
3468      return(p);
3469    }
3470 }
3471 // out of while when p==0 or p == const
3472 return(p);
3473}
3474
3475
3476
3477
3478//secundary procedures, all static
3479
3480static proc getlpCoeffs(poly q, poly p)
3481"
3482"
3483{list R; poly m; intvec cq,t,lv,rv,bla;
3484 int n = attrib(basering,"lV"); int d = attrib(basering,"uptodeg");
3485 int i;
3486 m = p/q;
3487 cq = leadexp(m);
3488 for (i = 1; i<= d; i++)
3489 {bla = cq[((i-1)*n+1)..(i*n)];
3490  if (bla == 0) {lv = cq[1..i*n]; cq = cq[(i*n+1)..(d*n)]; break;}
3491 }
3492
3493 d = size(cq) div n;
3494 for (i = 1; i<= d; i++)
3495 {bla = cq[((i-1)*n+1)..(i*n)];
3496  if (bla <> 0){rv = cq[((i-1)*n+1)..(d*n)]; break;}
3497 }
3498 return(list(monomial(lv),monomial(rv)));
3499}
3500
3501static proc lpReduce(poly p, poly g, int s)
3502"NOTE: shift can not exceed the degree bound, because s*g | p
3503"
3504{poly l,r,qt; int i;
3505 g = shiftPoly(g,s);
3506 list K = getlpCoeffs(lead(g),lead(p));
3507 l = K[1]; r = K[2];
3508 kill K;
3509 for (i = 1; i <= size(g); i++)
3510 {qt = qt + lpMult(lpMult(l,g[i]),r);
3511 }
3512 return((leadcoef(qt)*p - leadcoef(p)*qt));
3513}
3514
3515
3516static proc lpShrink(poly p)
3517"
3518"
3519{int n;
3520 if (typeof(attrib(basering,"isLetterplaceRing"))=="int")
3521 {n = attrib(basering,"lV");
3522  return(system("shrinktest",p,n));
3523 }
3524 else {ERROR("Basering is not a Letterplace ring!");}
3525}
3526
3527static proc checkAssumptions(poly p, ideal G)
3528"
3529"
3530{checkLPRing();
3531 checkAssumptionPoly(p);
3532 checkAssumptionIdeal(G);
3533 return();
3534}
3535
3536static proc checkLPRing();
3537"
3538"
3539{if (typeof(attrib(basering,"isLetterplaceRing"))=="string") {ERROR("Basering is not a Letterplace ring!");}
3540 return();
3541}
3542
3543static proc checkAssumptionIdeal(ideal G)
3544"PURPOSE:Check if all elements of ideal are elements of V'
3545"
3546{ideal L = lead(normalize(G));
3547 int i;
3548 for (i = 1; i <= ncols(G); i++) {if (!isContainedInVp(G[i])) {ERROR("Ideal containes elements not contained in V'");}}
3549 return();
3550}
3551
3552static proc checkAssumptionPoly(poly p)
3553"PURPOSE:Check if p is an element of V'
3554"
3555{poly l = lead(normalize(p));
3556 if (!isContainedInVp(l)) {ERROR("Polynomial is not contained in V'");}
3557 return();
3558}
3559
3560static proc isContainedInVp(poly p)
3561"PURPOSE: Check monomial for holes in the places
3562"
3563{int r = 0; intvec w;
3564 intvec l = leadexp(p);
3565 int n = attrib(basering,"lV"); int d = attrib(basering,"uptodeg");
3566 int i,j,c,c1;
3567 while (1 <= d)
3568 {w = l[1..n];
3569  if (w<>intvec(0)){break;}
3570   else {l = l[(n+1)..(n*d)]; d = d-1;}
3571 }
3572
3573 while (1 <= d)
3574  {for (j = 1; j <= n; j++)
3575   {if (l[j]<>0)
3576    {if (c1<>0){return(0);}
3577     if (c<>0){return(0);}
3578     if (l[j]<>1){return(0);}
3579     c=1;
3580    }
3581   }
3582   if (c == 0){c1=1;if (1 < d){l = l[(n+1)..(n*d)]; d = d-1;} else {d = d -1;}}
3583    else {c = 0; if (1 < d){l = l[(n+1)..(n*d)]; d = d-1;} else {d = d -1;}}
3584  }
3585 return(1);
3586}
3587
3588// under development for Roberto
3589static proc extractLinearPart(module M)
3590{
3591  /* returns vectors from a module whose max leadexp is 1 */
3592  /* does not take place nonlinearity into account yet */
3593  /* use rather kernel function isinV to get really nonlinear things */
3594  int i; int s = ncols(M);
3595  int answer = 1;
3596  vector v; module Ret;
3597  for(i=1; i<=s; i++)
3598  {
3599    if ( isLinearVector(M[i]) )
3600    {
3601      Ret = Ret, M[i];
3602    }
3603  }
3604  Ret = simplify(Ret,2);
3605  return(Ret);
3606}
3607
3608// under development for Roberto
3609static proc isLinearVector(vector v)
3610{
3611  /* returns true iff max leadexp is 1 */
3612  int i,j,k;
3613  intvec w;
3614  int s = size(v);
3615  poly p;
3616  int answer = 1;
3617  for(i=1; i<=s; i++)
3618  {
3619    p = v[i];
3620    while (p != 0)
3621    {
3622      w = leadexp(p);
3623      j = Max(w);
3624      if (j >=2)
3625      {
3626        answer = 0;
3627        return(answer);
3628      }
3629      p = p-lead(p);
3630    }
3631  }
3632  return(answer);
3633}
3634
3635
3636// // the following is to determine a shift of a mono/poly from the
3637// // interface
3638
3639// static proc whichshift(poly p, int numvars)
3640// {
3641// // numvars = number of vars of the orig free algebra
3642// // assume: we are in the letterplace ring
3643// // takes  monomial on the input
3644// poly q = lead(p);
3645// intvec v = leadexp(v);
3646// if (v==0) { return(int(0)); }
3647// int sv = size(v);
3648// int i=1;
3649// while ( (v[i]==0) && (i<sv) ) { i++; }
3650// i = sv div i;
3651// return(i);
3652// }
3653
3654
3655// LIB "qhmoduli.lib";
3656// static proc polyshift(poly p,  int numvars)
3657// {
3658//   poly q = p; int i = 0;
3659//   while (q!=0)
3660//   {
3661//     i = Max(i, whichshift(q,numvars));
3662//     q = q - lead(q);
3663//   }
3664//   return(q);
3665// }
3666
3667static proc lpAssumeViolation()
3668{
3669  // checks whether the global vars
3670  // uptodeg and lV are defined
3671  // returns Boolean : yes/no [for assume violation]
3672  def lpring = attrib(basering,"isLetterplaceRing");
3673  if ( typeof(lpring)!="int" )
3674  {
3675    //  if ( typeof(lpring)=="string" ) ??
3676    // basering is NOT lp Ring
3677    return(1);
3678  }
3679  def uptodeg = attrib(basering,"uptodeg");
3680  if ( typeof(uptodeg)!="int" )
3681  {
3682    return(1);
3683  }
3684  def lV = attrib(basering,"lV");
3685  if ( typeof(lV)!="int" )
3686  {
3687    return(1);
3688  }
3689  //  int i = ( defined(uptodeg) && (defined(lV)) );
3690  //  return ( !i );
3691  return(0);
3692}
3693
3694static proc bugSKing()
3695{
3696  LIB "freegb.lib";
3697  ring r=0,(a,b),dp;
3698  def R = makeLetterplaceRing(5);
3699  setring R;
3700  poly p = a(1);
3701  poly q = b(1);
3702  poly p2 = lpPower(p,2);
3703  lpMult(p2+q,q)-lpMult(p2,q)-lpMult(q,q); // now its 0
3704}
3705
3706static proc bugRucker()
3707{
3708  // needs unstatic lpMultX
3709  LIB "freegb.lib";
3710  ring r=0,(a,b,c,d,p,q,r,s,t,u,v,w),(a(7,1,1,7),dp);
3711  def R=makeLetterplaceRing(20,1);
3712  setring R;
3713  option(redSB); option(redTail);
3714  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);
3715  poly ttt = a(1)*v(2)*w(3)-p(1)*q(2)*r(3)*s(4)*t(5)*u(6)*d(7);
3716  // with lpMult
3717  lpMult(I[1],d(1)) - lpMult(a(1),I[2]); // spoly; has been incorrect before
3718  _ - ttt;
3719  // with lpMultX
3720  lpMultX(I[1],d(1)) - lpMultX(a(1),I[2]); // spoly; has been incorrect before
3721  _ - ttt;
3722}
3723
3724static proc checkWeightedExampleLP()
3725{
3726  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);
3727  def R = setLetterplaceAttributes(r,4,2); // supply R with letterplace structure
3728  setring R;
3729  poly a = x(1)*y(2)+x(1)+y(1); poly b = y(1)+3;
3730  lpMultX(b,a);
3731  lpMultX(a,b); // seems to work properly
3732}
3733
3734/* THE FOLLOWING ARE UNDER DEVELOPMENT
3735// copied following from freegb_wrkcp.lib by Karim Abou Zeid on 07.04.2017:
3736// makeLetterplaceRingElim(int d)
3737// makeLetterplaceRingNDO(int d)
3738// setLetterplaceAttributesElim(def R, int uptodeg, int lV)
3739// lpElimIdeal(ideal I)
3740// makeLetterplaceRingWt(int d, intvec W)
3741
3742static proc makeLetterplaceRingElim(int d)
3743"USAGE:  makeLetterplaceRingElim(d); d integers
3744RETURN:  ring
3745PURPOSE: creates a ring with an elimination ordering
3746NOTE: the matrix for the ordering looks as follows: first row is 1,..,0,1,0,..
3747@* then 0,1,0,...,0,0,1,0... and so on, lastly its lp
3748@* this ordering is only correct if only polys with same shift are compared
3749EXAMPLE: example makeLetterplaceRingElim; shows examples
3750"
3751{
3752
3753  // ToDo future: inherit positive weights in the orig ring
3754  // complain on nonpositive ones
3755
3756  // d = up to degree, will be shifted to d+1
3757  if (d<1) {"bad d"; return(0);}
3758
3759  int uptodeg = d; int lV = nvars(basering);
3760
3761  int ppl = printlevel-voice+2;
3762  string err = "";
3763
3764  int i,j,s; intvec iV,iVl;
3765  def save = basering;
3766  int D = d-1;
3767  list LR  = ringlist(save);
3768  list L, tmp, tmp2, tmp3;
3769  L[1] = LR[1]; // ground field
3770  L[4] = LR[4]; // quotient ideal
3771  tmp  = LR[2]; // varnames
3772  s = size(LR[2]);
3773  for (i=1; i<=D; i++)
3774  {
3775    for (j=1; j<=s; j++)
3776    {
3777      tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")";
3778    }
3779  }
3780  for (i=1; i<=s; i++)
3781  {
3782    tmp[i] = string(tmp[i])+"("+string(1)+")";
3783  }
3784  L[2] = tmp;
3785  L[3] = list();
3786  list OrigNames = LR[2];
3787  s = size(LR[3]);
3788  //creation of first block
3789
3790  if (s==2)
3791  {
3792    // not a blockord, 1 block + module ord
3793    tmp = LR[3][s]; // module ord
3794    for (i = 1; i <= lV;  i++)
3795    {
3796      iV = (0: lV);
3797      iV[i] = 1;
3798      iVl = iV;
3799      for (j = 1; j <= D; j++)
3800       { iVl = iVl,iV; }
3801      L[3][i] = list("a",iVl);
3802    }
3803//    for (i=1; i<=d; i++)
3804//    {
3805//      LR[3][s-1+i] = LR[3][1];
3806//    }
3807    //    LR[3][s+D] = tmp;
3808    //iV = (1:(d*lV));
3809    L[3][lV+1] = list("lp",(1:(d*lV)));
3810    L[3][lV+2] = tmp;
3811  }
3812  else {ERROR("Please set the ordering of basering to dp");}
3813//  if (s>2)
3814//  {
3815//    // there are s-1 blocks
3816//    int nb = s-1;
3817//    tmp = LR[3][s]; // module ord to place at the very end
3818//   tmp2 = LR[3]; tmp2 = tmp2[1..nb];
3819//    LR[3][1] = list("a",LTO);
3820//    //tmp3[1] = list("a",intvec(1: int(d*lV))); // deg-ord, insert as the 1st
3821//    for (i=1; i<=d; i++)
3822//    {
3823//      tmp3 = tmp3 + tmp2;
3824//    }
3825//    tmp3 = tmp3 + list(tmp);
3826//    LR[3] = tmp3;
3827//     for (i=1; i<=d; i++)
3828//     {
3829//       for (j=1; j<=nb; j++)
3830//       {
3831//         //        LR[3][i*nb+j+1]= LR[3][j];
3832//         LR[3][i*nb+j+1]= tmp2[j];
3833//       }
3834//     }
3835//     //    size(LR[3]);
3836//     LR[3][(s-1)*d+2] = tmp;
3837//     LR[3] = list("a",intvec(1: int(d*lV))) + LR[3]; // deg-ord, insert as the 1st
3838    // remove everything behind nb*(D+1)+1 ?
3839    //    tmp = LR[3];
3840    //    LR[3] = tmp[1..size(tmp)-1];
3841 // }
3842 // L[3] = LR[3];
3843  def @R = ring(L);
3844  //  setring @R;
3845  //  int uptodeg = d; int lV = nvars(basering); // were defined before
3846  def @@R = setLetterplaceAttributesElim(@R,uptodeg,lV);
3847  return (@@R);
3848}
3849example
3850{
3851  "EXAMPLE:"; echo = 2;
3852  ring r = 0,(x,y,z),lp;
3853  def A = makeLetterplaceRingElim(2);
3854  setring A;
3855  A;
3856  attrib(A,"isLetterplaceRing");
3857  attrib(A,"uptodeg");  // degree bound
3858  attrib(A,"lV"); // number of variables in the main block
3859}
3860
3861
3862
3863static proc makeLetterplaceRingNDO(int d)
3864"USAGE:  makeLetterplaceRingNDO(d); d an integer
3865RETURN:  ring
3866PURPOSE: creates a ring with a non-degree first ordering, suitable for
3867@* the use of non-homogeneous letterplace
3868NOTE: the matrix for the ordering looks as follows:
3869@*    'd' blocks of shifted original variables
3870EXAMPLE: example makeLetterplaceRingNDO; shows examples
3871"
3872{
3873
3874  // ToDo future: inherit positive weights in the orig ring
3875  // complain on nonpositive ones
3876
3877  // d = up to degree, will be shifted to d+1
3878  if (d<1) {"bad d"; return(0);}
3879
3880  int uptodeg = d; int lV = nvars(basering);
3881
3882  int ppl = printlevel-voice+2;
3883  string err = "";
3884
3885  int i,j,s;
3886  def save = basering;
3887  int D = d-1;
3888  list LR  = ringlist(save);
3889  list L, tmp, tmp2, tmp3;
3890  L[1] = LR[1]; // ground field
3891  L[4] = LR[4]; // quotient ideal
3892  tmp  = LR[2]; // varnames
3893  s = size(LR[2]);
3894  for (i=1; i<=D; i++)
3895  {
3896    for (j=1; j<=s; j++)
3897    {
3898      tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")";
3899    }
3900  }
3901  for (i=1; i<=s; i++)
3902  {
3903    tmp[i] = string(tmp[i])+"("+string(1)+")";
3904  }
3905  L[2] = tmp;
3906  list OrigNames = LR[2];
3907  // ordering: one 1..1 a above
3908  // ordering: d blocks of the ord on r
3909  // try to get whether the ord on r is blockord itself
3910  // TODO: make L(2) ordering! exponent is maximally 2
3911  s = size(LR[3]);
3912  if (s==2)
3913  {
3914    // not a blockord, 1 block + module ord
3915    tmp = LR[3][s]; // module ord
3916    for (i=1; i<=d; i++)
3917    {
3918      LR[3][i] = LR[3][1];
3919    }
3920    //    LR[3][s+D] = tmp;
3921    LR[3][d+1] = tmp;
3922    //LR[3][1] = list("a",intvec(1: int(d*lV))); // deg-ord not wanted here
3923  }
3924  if (s>2)
3925  {
3926    // there are s-1 blocks
3927    int nb = s-1;
3928    tmp = LR[3][s]; // module ord to place at the very end
3929    tmp2 = LR[3]; tmp2 = tmp2[1..nb];
3930    //tmp3[1] = list("a",intvec(1: int(d*lV))); // deg-ord not wanted here
3931    for (i=1; i<=d; i++)
3932    {
3933      tmp3 = tmp3 + tmp2;
3934    }
3935    tmp3 = tmp3 + list(tmp);
3936    LR[3] = tmp3;
3937//     for (i=1; i<=d; i++)
3938//     {
3939//       for (j=1; j<=nb; j++)
3940//       {
3941//         //        LR[3][i*nb+j+1]= LR[3][j];
3942//         LR[3][i*nb+j+1]= tmp2[j];
3943//       }
3944//     }
3945//     //    size(LR[3]);
3946//     LR[3][(s-1)*d+2] = tmp;
3947//     LR[3] = list("a",intvec(1: int(d*lV))) + LR[3]; // deg-ord, insert as the 1st
3948    // remove everything behind nb*(D+1)+1 ?
3949    //    tmp = LR[3];
3950    //    LR[3] = tmp[1..size(tmp)-1];
3951  }
3952  L[3] = LR[3];
3953  def @R = ring(L);
3954  //  setring @R;
3955  //  int uptodeg = d; int lV = nvars(basering); // were defined before
3956  def @@R = setLetterplaceAttributes(@R,uptodeg,lV);
3957  return (@@R);
3958}
3959example
3960{
3961  "EXAMPLE:"; echo = 2;
3962  ring r = 0,(x,y,z),lp;
3963  def A = makeLetterplaceRingNDO(2);
3964  setring A;
3965  A;
3966  attrib(A,"isLetterplaceRing");
3967  attrib(A,"uptodeg");  // degree bound
3968  attrib(A,"lV"); // number of variables in the main block
3969}
3970
3971static proc setLetterplaceAttributesElim(def R, int uptodeg, int lV)
3972"USAGE: setLetterplaceAttributesElim(R, d, b, eV); R a ring, b,d, eV integers
3973RETURN: ring with special attributes set
3974PURPOSE: sets attributes for a letterplace ring:
3975@*      'isLetterplaceRing' = true, 'uptodeg' = d, 'lV' = b, 'eV' = eV, where
3976@*      'uptodeg' stands for the degree bound,
3977@*      'lV' for the number of variables in the block 0
3978@*      'eV' for the number of elimination variables
3979NOTE: Activate the resulting ring by using @code{setring}
3980"
3981{
3982  if (uptodeg*lV != nvars(R))
3983  {
3984    ERROR("uptodeg and lV do not agree on the basering!");
3985  }
3986
3987
3988    // Set letterplace-specific attributes for the output ring!
3989  attrib(R, "uptodeg", uptodeg);
3990  attrib(R, "lV", lV);
3991  attrib(R, "isLetterplaceRing", 1);
3992  attrib(R, "HasElimOrd", 1);
3993  return (R);
3994}
3995example
3996{
3997  "EXAMPLE:"; echo = 2;
3998  ring r = 0,(x(1),y(1),x(2),y(2),x(3),y(3),x(4),y(4)),dp;
3999  def R = setLetterplaceAttributesElim(r, 4, 2, 1); setring R;
4000  attrib(R,"isLetterplaceRing");
4001  lieBracket(x(1),y(1),2);
4002}
4003
4004
4005static proc lpElimIdeal(ideal I)
4006"
4007does not work for degree reasons (deg function does not work for lp rings -> newone!)
4008"
4009{
4010  def lpring = attrib(basering,"isLetterplaceRing");
4011  def lpEO =  attrib(basering,"HasElimOrd");
4012  if ( typeof(lpring)!="int" && typeof(lpEO)!="int")
4013  {
4014    ERROR("Ring is not a lp-ring with an elimination ordering");
4015  }
4016
4017  //int nE = attrib(basering, "eV");
4018
4019  return(letplaceGBasis(I));
4020}
4021
4022
4023static proc makeLetterplaceRingWt(int d, intvec W)
4024"USAGE:  makeLetterplaceRingWt(d,W); d an integer, W a vector of positive integers
4025RETURN:  ring
4026PURPOSE: creates a ring with a special ordering, suitable for
4027@* the use of non-homogeneous letterplace
4028NOTE: the matrix for the ordering looks as follows: first row is W,W,W,...
4029@* then there come 'd' blocks of shifted original variables
4030EXAMPLE: example makeLetterplaceRing2; shows examples
4031"
4032{
4033
4034  // ToDo future: inherit positive weights in the orig ring
4035  // complain on nonpositive ones
4036
4037  // d = up to degree, will be shifted to d+1
4038  if (d<1) {"bad d"; return(0);}
4039
4040  int uptodeg = d; int lV = nvars(basering);
4041
4042  //check weightvector
4043  if (size(W) <> lV) {"bad weights"; return(0);}
4044
4045  int i;
4046  for (i = 1; i <= size(W); i++) {if (W[i] < 0) {"bad weights"; return(0);}}
4047  intvec Wt = W;
4048  for (i = 2; i <= d; i++) {Wt = Wt, W;}
4049  kill i;
4050
4051  int ppl = printlevel-voice+2;
4052  string err = "";
4053
4054  int i,j,s;
4055  def save = basering;
4056  int D = d-1;
4057  list LR  = ringlist(save);
4058  list L, tmp, tmp2, tmp3;
4059  L[1] = LR[1]; // ground field
4060  L[4] = LR[4]; // quotient ideal
4061  tmp  = LR[2]; // varnames
4062  s = size(LR[2]);
4063  for (i=1; i<=D; i++)
4064  {
4065    for (j=1; j<=s; j++)
4066    {
4067      tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")";
4068    }
4069  }
4070  for (i=1; i<=s; i++)
4071  {
4072    tmp[i] = string(tmp[i])+"("+string(1)+")";
4073  }
4074  L[2] = tmp;
4075  list OrigNames = LR[2];
4076  // ordering: one 1..1 a above
4077  // ordering: d blocks of the ord on r
4078  // try to get whether the ord on r is blockord itself
4079  // TODO: make L(2) ordering! exponent is maximally 2
4080  s = size(LR[3]);
4081  if (s==2)
4082  {
4083    // not a blockord, 1 block + module ord
4084    tmp = LR[3][s]; // module ord
4085    for (i=1; i<=d; i++)
4086    {
4087      LR[3][s-1+i] = LR[3][1];
4088    }
4089    //    LR[3][s+D] = tmp;
4090    LR[3][s+1+D] = tmp;
4091    LR[3][1] = list("a",Wt); // deg-ord
4092  }
4093  if (s>2)
4094  {
4095    // there are s-1 blocks
4096    int nb = s-1;
4097    tmp = LR[3][s]; // module ord to place at the very end
4098    tmp2 = LR[3]; tmp2 = tmp2[1..nb];
4099    tmp3[1] = list("a",Wt); // deg-ord, insert as the 1st
4100    for (i=1; i<=d; i++)
4101    {
4102      tmp3 = tmp3 + tmp2;
4103    }
4104    tmp3 = tmp3 + list(tmp);
4105    LR[3] = tmp3;
4106
4107  }
4108  L[3] = LR[3];
4109  def @R = ring(L);
4110  //  setring @R;
4111  //  int uptodeg = d; int lV = nvars(basering); // were defined before
4112  def @@R = setLetterplaceAttributes(@R,uptodeg,lV);
4113  return (@@R);
4114}
4115example
4116{
4117  "EXAMPLE:"; echo = 2;
4118  ring r = 0,(x,y,z),(dp(1),dp(2));
4119  def A = makeLetterplaceRingWt(2,intvec(1,2,3));
4120  setring A;
4121  A;
4122  attrib(A,"isLetterplaceRing");
4123  attrib(A,"uptodeg");  // degree bound
4124  attrib(A,"lV"); // number of variables in the main block
4125}
4126*/
Note: See TracBrowser for help on using the repository browser.