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

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