source:git/Singular/LIB/freegb.lib@320db1

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