source: git/Singular/LIB/freegb.lib @ e7d857

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