source: git/Singular/LIB/freegb.lib @ 291c20

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