# source:git/Singular/LIB/freegb.lib@db0c264

spielwiese
Last change on this file since db0c264 was db0c264, checked in by Viktor Levandovskyy <levandov@…>, 15 years ago
*levandov: letterplace changes related to docu git-svn-id: file:///usr/local/Singular/svn/trunk@11272 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 48.9 KB
Line
1//////////////////////////////////////////////////////////////////////////////
2version="\$Id: freegb.lib,v 1.13 2008-12-26 18:15:40 levandov Exp \$";
3category="Noncommutative";
4info="
5LIBRARY: freegb.lib   Twosided Noncommutative Groebner bases in Free Algebras
6AUTHOR: Viktor Levandovskyy,     levandov@math.rwth-aachen.de
7
8PROCEDURES:
9freegbRing(d);    creates a ring with d blocks of shifted original variables
10freegbasis(L, int n);   compute two-sided Groebner basis of ideal, encoded via L, up to degree n
11
12AUXILIARY PROCEDURES:
13
14lp2lstr(K, s);      convert letter-place ideal to a list of modules
15lst2str(L[, n]);     convert a list (of modules) into polynomials in free algebra
16mod2str(M[, n]); convert a module into a polynomial in free algebra
17vct2str(M[, n]);   convert a vector into a word in free algebra
18Liebr(a,b[, N]);    compute Lie bracket ab-ba of two letterplace polynomials
19Serre(A,z);          compute the ideal of Serre's relations associated to a generalized Cartan matrix A
20isVar(p);             check whether p is a power of a single variable
21
23"
24
25// this library computes two-sided GB of an ideal
26// in a free associative algebra
27
28// a monomial is encoded via a vector V
29// where V[1] = coefficient
30// V[1+i] = the corresponding symbol
31
32LIB "discretize.lib"; // for replace
33LIB "qhmoduli.lib"; // for Max
34
35proc testfreegblib()
36{
37  example freegbRing;
38  example freegbasis;
39    "AUXILIARY PROCEDURES: ";
40  example lp2lstr;
41  example lst2str;
42  example mod2str;
43  example vct2str;
44  example Liebr;
45  example Serre;
46  example isVar;
47}
48
49
50// obsolete?
51
52static proc lshift(module M, int s, string varing, def lpring)
53{
54  // FINALLY IMPLEMENTED AS A PART OT THE CODE
55  // shifts a poly from the ring @R to s positions
56  // M lives in varing, the result in lpring
57  // to be run from varing
58  int i, j, k, sm, sv;
59  vector v;
60  //  execute("setring "+lpring);
61  setring lpring;
62  poly @@p;
63  ideal I;
64  execute("setring "+varing);
65  sm = ncols(M);
66  for (i=1; i<=s; i++)
67  {
68    // modules, e.g. free polynomials
69    for (j=1; j<=sm; j++)
70    {
71      //vectors, e.g. free monomials
72      v  = M[j];
73      sv = size(v);
74      sp = "@@p = @@p + ";
75      for (k=2; k<=sv; k++)
76      {
77        sp = sp + string(v[k])+"("+string(k-1+s)+")*";
78      }
79      sp = sp + string(v[1])+";"; // coef;
80      setring lpring;
81      //      execute("setring "+lpring);
82      execute(sp);
83      execute("setring "+varing);
84    }
85    setring lpring;
86    //    execute("setring "+lpring);
87    I = I,@@p;
88    @@p = 0;
89  }
90  setring lpring;
91  //execute("setring "+lpring);
92  export(I);
93  //  setring varing;
94  execute("setring "+varing);
95}
96
97static proc skip0(vector v)
98{
99  // skips zeros in a vector, producing another vector
100  int sv = nrows(v);
101  int sw = size(v);
102  if (sv == sw)
103  {
104    return(v);
105  }
106  int i;
107  int j=1;
108  vector w;
109  for (i=1; i<=sv; i++)
110  {
111    if (v[i] != 0)
112    {
113      w = w + v[i]*gen(j);
114      j++;
115    }
116  }
117  return(w);
118}
119
120proc lst2str(list L, list #)
121"USAGE:  lst2str(L[,n]);  L a list of modules, n an optional integer
122RETURN:  list (of strings)
123PURPOSE: convert a list (of modules) into polynomials in free algebra
124EXAMPLE: example lst2str; shows examples
125NOTE: if an optional integer is not 0, stars signs are used in multiplication
126"
127{
128  // returns a list of strings
129  // being sentences in words built from L
130  // if #[1] = 1, use * between generators
131  int useStar = 0;
132  if ( size(#)>0 )
133  {
134    if (#[1])
135    {
136      useStar = 1;
137    }
138  }
139  int i;
140  int s    = size(L);
141  list N;
142  for(i=1; i<=s; i++)
143  {
144    if ((typeof(L[i]) == "module") || (typeof(L[i]) == "matrix") )
145    {
146      N[i] = mod2str(L[i],useStar);
147    }
148    else
149    {
150      "module or matrix expected in the list";
151      return(N);
152    }
153  }
154  return(N);
155}
156example
157{
158  "EXAMPLE:"; echo = 2;
159  ring r = 0,(x,y,z),(dp(1),dp(2));
160  module M = [-1,x,y],[-7,y,y],[3,x,x];
161  module N = [1,x,y,x,y],[-2,y,x,y,x],[6,x,y,y,x,y];
162  list L; L[1] = M; L[2] = N;
163  lst2str(L);
164  lst2str(L[1],1);
165}
166
167
168proc mod2str(module M, list #)
169"USAGE:  mod2str(M[,n]);  M a module, n an optional integer
170RETURN:  string
171PURPOSE: convert a module into a polynomial in free algebra
172EXAMPLE: example mod2str; shows examples
173NOTE: if an optional integer is not 0, stars signs are used in multiplication
174"
175{
176  // returns a string
177  // a sentence in words built from M
178  // if #[1] = 1, use * between generators
179  int useStar = 0;
180  if ( size(#)>0 )
181  {
182    if (#[1])
183    {
184      useStar = 1;
185    }
186  }
187  int i;
188  int s    = ncols(M);
189  string t;
190  string mp;
191  for(i=1; i<=s; i++)
192  {
193    mp = vct2str(M[i],useStar);
194    if (mp[1] == "-")
195    {
196      t = t + mp;
197    }
198    else
199    {
200      t = t + "+" + mp;
201    }
202  }
203  if (t[1]=="+")
204  {
205    t = t[2..size(t)]; // remove first "+"
206  }
207  return(t);
208}
209example
210{
211  "EXAMPLE:"; echo = 2;
212  ring r = 0,(x,y,z),(dp);
213  module M = [1,x,y,x,y],[-2,y,x,y,x],[6,x,y,y,x,y];
214  mod2str(M);
215  mod2str(M,1);
216}
217
218proc vct2str(vector v, list #)
219"USAGE:  vct2str(v[,n]);  v a vector, n an optional integer
220RETURN:  string
221PURPOSE: convert a vector into a word in free algebra
222EXAMPLE: example vct2str; shows examples
223NOTE: if an optional integer is not 0, stars signs are used in multiplication
224"
225{
226  // if #[1] = 1, use * between generators
227  int useStar = 0;
228  if ( size(#)>0 )
229  {
230    if (#[1])
231    {
232      useStar = 1;
233    }
234  }
235  int ppl = printlevel-voice+2;
236  // for a word, encoded by v
237  // produces a string for it
238  v = skip0(v);
240  int s = size(v);
241  string vs,vv,vp,err;
242  int i,j,p,q;
243  for (i=1; i<=s-1; i++)
244  {
245    p     = isVar(v[i+1]);
246    if (p==0)
247    {
248      err = "Error: monomial expected at" + string(i+1);
249      dbprint(ppl,err);
250      return("_");
251    }
252    if (p==1)
253    {
254      if (useStar && (size(vs) >0))       {   vs = vs + "*"; }
255        vs = vs + string(v[i+1]);
256    }
257    else //power
258    {
259      vv = string(v[i+1]);
260      q = find(vv,"^");
261      if (q==0)
262      {
263        q = find(vv,string(p));
264        if (q==0)
265        {
266          err = "error in find for string "+vv;
267          dbprint(ppl,err);
268          return("_");
269        }
270      }
271      // q>0
272      vp = vv[1..q-1];
273      for(j=1;j<=p;j++)
274      {
275         if (useStar && (size(vs) >0))       {   vs = vs + "*"; }
276         vs = vs + vp;
277      }
278    }
279  }
280  string scf;
281  if (cf == -1)
282  {
283    scf = "-";
284  }
285  else
286  {
287    scf = string(cf);
288    if (cf == 1)
289    {
290      scf = "";
291    }
292  }
293  if (useStar && (size(scf) >0) && (scf!="-") )       {   scf = scf + "*"; }
294  vs = scf + vs;
295  return(vs);
296}
297example
298{
299  "EXAMPLE:"; echo = 2;
300  ring r = (0,a),(x,y3,z(1)),dp;
301  vector v = [-7,x,y3^4,x2,z(1)^3];
302  vct2str(v);
303  vct2str(v,1);
304  vector w = [-7a^5+6a,x,y3,y3,x,z(1),z(1)];
305  vct2str(w);
306  vct2str(w,1);
307}
308
309proc isVar(poly p)
310"USAGE:  isVar(p);  poly p
311RETURN:  int
312PURPOSE: checks whether p is a power of a single variable from the basering.
313@* Returns the exponent or 0 is p is not a power of  a single variable.
314EXAMPLE: example isVar; shows examples
315"
316{
317  // checks whether p is a variable indeed
318  // if it's a power of a variable, returns the power
319  if (p==0) {  return(0); } //"p=0";
321  if ( (p-lead(p)) !=0 ) { return(0); } // "p-lm(p)>0";
323  int s = size(v);
324  int i=1;
325  int cnt = 0;
326  int pwr = 0;
327  for (i=1; i<=s; i++)
328  {
329    if (v[i] != 0)
330    {
331      cnt++;
332      pwr = v[i];
333    }
334  }
335  //  "cnt:";  cnt;
336  if (cnt==1) { return(pwr); }
337  else { return(0); }
338}
339example
340{
341  "EXAMPLE:"; echo = 2;
342  ring r = 0,(x,y),dp;
343  poly f = xy+1;
344  isVar(f);
345  poly g = xy;
346  isVar(g);
347  poly h = y^3;
348  isVar(h);
349  poly i = 1;
350  isVar(i);
351}
352
353// new conversion routines
354
355static proc id2words(ideal I, int d)
356{
357  // NOT FINISHED
358  // input: ideal I of polys in letter-place notation
359  // in the ring with d real vars
360  // output: the list of strings: associative words
361  // extract names of vars
362  int i,m,n; string s; string place = "(1)";
363  list lv;
364  for(i=1; i<=d; i++)
365  {
366    s = string(var(i));
367    // get rid of place
368    n = find(s, place);
369    if (n>0)
370    {
371      s = s[1..n-1];
372    }
373    lv[i] = s;
374  }
375  poly p,q;
376  for (i=1; i<=ncols(I); i++)
377  {
378    if (I[i] != 0)
379    {
380      p = I[i];
381      while (p!=0)
382      {
384      }
385    }
386  }
387
388  return(lv);
389}
390example
391{
392  "EXAMPLE:"; echo = 2;
393  ring r = 0,(x(1),y(1),z(1),x(2),y(2),z(2)),dp;
394  ideal I = x(1)*y(2) -z(1)*x(2);
395  id2words(I,3);
396}
397
398static proc mono2word(poly p, int d)
399{
400}
401
402// given the element -7xy^2x, it is represented as [-7,x,y^2,x] or as [-7,x,y,y,x]
403// use the orig ord on (x,y,z) and expand it blockwise to (x(i),y(i),z(i))
404
405// the correspondences:
406// monomial in K<x,y,z>    <<--->> vector in R
407// polynomial in K<x,y,z>  <<--->> list of vectors (matrix/module) in R
408// ideal in K<x,y,z>       <<--->> list of matrices/modules in R
409
410
411// 1. form a new ring
412// 2. NOP
413// 3. compute GB -> with the kernel stuff
414// 4. skip shifted elts (check that no such exist?)
415// 5. go back to orig vars, produce strings/modules
416// 6. return the result
417
418proc freegbasis(list LM, int d)
419"USAGE:  freegbasis(L, d);  L a list of modules, d an integer
420RETURN:  ring
421PURPOSE: compute the two-sided Groebner basis of an ideal, encoded by L in
422the free associative algebra, up to degree d
423EXAMPLE: example freegbasis; shows examples
424"
425{
426  // d = up to degree, will be shifted to d+1
427  if (d<1) {"bad d"; return(0);}
428
429  int ppl = printlevel-voice+2;
430  string err = "";
431
432  int i,j,s;
433  def save = basering;
434  // determine max no of places in the input
435  int slm = size(LM); // numbers of polys in the ideal
436  int sm;
437  intvec iv;
438  module M;
439  for (i=1; i<=slm; i++)
440  {
441    // modules, e.g. free polynomials
442    M  = LM[i];
443    sm = ncols(M);
444    for (j=1; j<=sm; j++)
445    {
446      //vectors, e.g. free monomials
447      iv = iv, size(M[j])-1; // 1 place is reserved by the coeff
448    }
449  }
450  int D = Max(iv); // max size of input words
451  if (d<D) {"bad d"; return(LM);}
452  D = D + d-1;
453  //  D = d;
454  list LR  = ringlist(save);
455  list L, tmp;
456  L[1] = LR[1]; // ground field
457  L[4] = LR[4]; // quotient ideal
458  tmp  = LR[2]; // varnames
459  s = size(LR[2]);
460  for (i=1; i<=D; i++)
461  {
462    for (j=1; j<=s; j++)
463    {
464      tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")";
465    }
466  }
467  for (i=1; i<=s; i++)
468  {
469    tmp[i] = string(tmp[i])+"("+string(1)+")";
470  }
471  L[2] = tmp;
472  list OrigNames = LR[2];
473  // ordering: d blocks of the ord on r
474  // try to get whether the ord on r is blockord itself
475  s = size(LR[3]);
476  if (s==2)
477  {
478    // not a blockord, 1 block + module ord
479    tmp = LR[3][s]; // module ord
480    for (i=1; i<=D; i++)
481    {
482      LR[3][s-1+i] = LR[3][1];
483    }
484    LR[3][s+D] = tmp;
485  }
486  if (s>2)
487  {
488    // there are s-1 blocks
489    int nb = s-1;
490    tmp = LR[3][s]; // module ord
491    for (i=1; i<=D; i++)
492    {
493      for (j=1; j<=nb; j++)
494      {
495        LR[3][i*nb+j] = LR[3][j];
496      }
497    }
498    //    size(LR[3]);
499    LR[3][nb*(D+1)+1] = tmp;
500  }
501  L[3] = LR[3];
502  def @R = ring(L);
503  setring @R;
504  ideal I;
505  poly @p;
506  s = size(OrigNames);
507  //  "s:";s;
508  // convert LM to canonical vectors (no powers)
509  setring save;
510  kill M; // M was defined earlier
511  module M;
512  slm = size(LM); // numbers of polys in the ideal
513  int sv,k,l;
514  vector v;
515  //  poly p;
516  string sp;
517  setring @R;
518  poly @@p=0;
519  setring save;
520  for (l=1; l<=slm; l++)
521  {
522    // modules, e.g. free polynomials
523    M  = LM[l];
524    sm = ncols(M); // in intvec iv the sizes are stored
525    // modules, e.g. free polynomials
526    for (j=1; j<=sm; j++)
527    {
528      //vectors, e.g. free monomials
529      v  = M[j];
530      sv = size(v);
531      //        "sv:";sv;
532      sp = "@@p = @@p + ";
533      for (k=2; k<=sv; k++)
534      {
535        sp = sp + string(v[k])+"("+string(k-1)+")*";
536      }
537      sp = sp + string(v[1])+";"; // coef;
538      setring @R;
539      execute(sp);
540      setring save;
541    }
542    setring @R;
543    //      "@@p:"; @@p;
544    I = I,@@p;
545    @@p = 0;
546    setring save;
547  }
548  kill sp;
549  // 3. compute GB
550  setring @R;
551  dbprint(ppl,"computing GB");
552  ideal J = system("freegb",I,d,nvars(save));
553  //  ideal J = slimgb(I);
554  dbprint(ppl,J);
555  // 4. skip shifted elts
556  ideal K = select1(J,1..s); // s = size(OrigNames)
557  dbprint(ppl,K);
558  dbprint(ppl, "done with GB");
559  // K contains vars x(1),...z(1) = images of originals
560  // 5. go back to orig vars, produce strings/modules
561  if (K[1] == 0)
562  {
563    "no reasonable output, GB gives 0";
564    return(0);
565  }
566  int sk = size(K);
567  int sp, sx, a, b;
568  intvec x;
569  poly p,q;
570  poly pn;
571  // vars in 'save'
572  setring save;
573  module N;
574  list LN;
575  vector V;
576  poly pn;
577  // test and skip exponents >=2
578  setring @R;
579  for(i=1; i<=sk; i++)
580  {
581    p  = K[i];
582    while (p!=0)
583    {
585      //      "processing q:";q;
587      sx = size(x);
588      for(k=1; k<=sx; k++)
589      {
590        if ( x[k] >= 2 )
591        {
592          err = "skip: the value x[k] is " + string(x[k]);
593          dbprint(ppl,err);
594          //            return(0);
595          K[i] = 0;
596          p    = 0;
597          q    = 0;
598          break;
599        }
600      }
601      p  = p - q;
602    }
603  }
604  K  = simplify(K,2);
605  sk = size(K);
606  for(i=1; i<=sk; i++)
607  {
608    //    setring save;
609    //    V  = 0;
610    setring @R;
611    p  = K[i];
612    while (p!=0)
613    {
615      err =  "processing q:" + string(q);
616      dbprint(ppl,err);
618      sx = size(x);
620      setring save;
621      pn = imap(@R,pn);
622      V  = V + leadcoef(pn)*gen(1);
623      for(k=1; k<=sx; k++)
624      {
625        if (x[k] ==1)
626        {
627          a = k / s; // block number=a+1, a!=0
628          b = k % s; // remainder
629          //          printf("a: %s, b: %s",a,b);
630          if (b == 0)
631          {
632            // that is it's the last var in the block
633            b = s;
634            a = a-1;
635          }
636          V = V + var(b)*gen(a+2);
637        }
638//         else
639//         {
640//           printf("error: the value x[k] is %s", x[k]);
641//           return(0);
642//         }
643      }
644      err = "V: " + string(V);
645      dbprint(ppl,err);
646      //      printf("V: %s", string(V));
647      N = N,V;
648      V  = 0;
649      setring @R;
650      p  = p - q;
651      pn = 0;
652    }
653    setring save;
654    LN[i] = simplify(N,2);
655    N     = 0;
656  }
657  setring save;
658  return(LN);
659}
660example
661{
662  "EXAMPLE:"; echo = 2;
663  ring r = 0,(x,y,z),(dp(1),dp(2));
664  module M = [-1,x,y],[-7,y,y],[3,x,x];
665  module N = [1,x,y,x],[-1,y,x,y];
666  list L; L[1] = M; L[2] = N;
667  lst2str(L);
668  def U = freegbasis(L,5);
669  lst2str(U);
670}
671
672static proc crs(list LM, int d)
673"USAGE:  crs(L, d);  L a list of modules, d an integer
674RETURN:  ring
675PURPOSE: create a ring and shift the ideal
676EXAMPLE: example crs; shows examples
677"
678{
679  // d = up to degree, will be shifted to d+1
680  if (d<1) {"bad d"; return(0);}
681
682  int ppl = printlevel-voice+2;
683  string err = "";
684
685  int i,j,s;
686  def save = basering;
687  // determine max no of places in the input
688  int slm = size(LM); // numbers of polys in the ideal
689  int sm;
690  intvec iv;
691  module M;
692  for (i=1; i<=slm; i++)
693  {
694    // modules, e.g. free polynomials
695    M  = LM[i];
696    sm = ncols(M);
697    for (j=1; j<=sm; j++)
698    {
699      //vectors, e.g. free monomials
700      iv = iv, size(M[j])-1; // 1 place is reserved by the coeff
701    }
702  }
703  int D = Max(iv); // max size of input words
704  if (d<D) {"bad d"; return(LM);}
705  D = D + d-1;
706  //  D = d;
707  list LR  = ringlist(save);
708  list L, tmp;
709  L[1] = LR[1]; // ground field
710  L[4] = LR[4]; // quotient ideal
711  tmp  = LR[2]; // varnames
712  s = size(LR[2]);
713  for (i=1; i<=D; i++)
714  {
715    for (j=1; j<=s; j++)
716    {
717      tmp[i*s+j] = string(tmp[j])+"("+string(i)+")";
718    }
719  }
720  for (i=1; i<=s; i++)
721  {
722    tmp[i] = string(tmp[i])+"("+string(0)+")";
723  }
724  L[2] = tmp;
725  list OrigNames = LR[2];
726  // ordering: d blocks of the ord on r
727  // try to get whether the ord on r is blockord itself
728  s = size(LR[3]);
729  if (s==2)
730  {
731    // not a blockord, 1 block + module ord
732    tmp = LR[3][s]; // module ord
733    for (i=1; i<=D; i++)
734    {
735      LR[3][s-1+i] = LR[3][1];
736    }
737    LR[3][s+D] = tmp;
738  }
739  if (s>2)
740  {
741    // there are s-1 blocks
742    int nb = s-1;
743    tmp = LR[3][s]; // module ord
744    for (i=1; i<=D; i++)
745    {
746      for (j=1; j<=nb; j++)
747      {
748        LR[3][i*nb+j] = LR[3][j];
749      }
750    }
751    //    size(LR[3]);
752    LR[3][nb*(D+1)+1] = tmp;
753  }
754  L[3] = LR[3];
755  def @R = ring(L);
756  setring @R;
757  ideal I;
758  poly @p;
759  s = size(OrigNames);
760  //  "s:";s;
761  // convert LM to canonical vectors (no powers)
762  setring save;
763  kill M; // M was defined earlier
764  module M;
765  slm = size(LM); // numbers of polys in the ideal
766  int sv,k,l;
767  vector v;
768  //  poly p;
769  string sp;
770  setring @R;
771  poly @@p=0;
772  setring save;
773  for (l=1; l<=slm; l++)
774  {
775    // modules, e.g. free polynomials
776    M  = LM[l];
777    sm = ncols(M); // in intvec iv the sizes are stored
778    for (i=0; i<=d-iv[l]; i++)
779    {
780      // modules, e.g. free polynomials
781      for (j=1; j<=sm; j++)
782      {
783        //vectors, e.g. free monomials
784        v  = M[j];
785        sv = size(v);
786        //        "sv:";sv;
787        sp = "@@p = @@p + ";
788        for (k=2; k<=sv; k++)
789        {
790          sp = sp + string(v[k])+"("+string(k-2+i)+")*";
791        }
792        sp = sp + string(v[1])+";"; // coef;
793        setring @R;
794        execute(sp);
795        setring save;
796      }
797      setring @R;
798      //      "@@p:"; @@p;
799      I = I,@@p;
800      @@p = 0;
801      setring save;
802    }
803  }
804  setring @R;
805  export I;
806  return(@R);
807}
808example
809{
810  "EXAMPLE:"; echo = 2;
811  ring r = 0,(x,y,z),(dp(1),dp(2));
812  module M = [-1,x,y],[-7,y,y],[3,x,x];
813  module N = [1,x,y,x],[-1,y,x,y];
814  list L; L[1] = M; L[2] = N;
815  lst2str(L);
816  def U = crs(L,5);
817  setring U; U;
818  I;
819}
820
821static proc polylen(ideal I)
822{
823  // returns the ideal of length of polys
824  int i;
825  intvec J;
826  number s = 0;
827  for(i=1;i<=ncols(I);i++)
828  {
829    J[i] = size(I[i]);
830    s = s + J[i];
831  }
832  printf("the sum of length %s",s);
833  //  print(s);
834  return(J);
835}
836
837proc freegbRing(int d)
838"USAGE:  freegbRing(d); d an integer
839RETURN:  ring
840PURPOSE: creates a ring with d blocks of shifted original variables
841EXAMPLE: example freegbRing; shows examples
842"
843{
844  // d = up to degree, will be shifted to d+1
845  if (d<1) {"bad d"; return(0);}
846
847  int ppl = printlevel-voice+2;
848  string err = "";
849
850  int i,j,s;
851  def save = basering;
852  int D = d-1;
853  list LR  = ringlist(save);
854  list L, tmp;
855  L[1] = LR[1]; // ground field
856  L[4] = LR[4]; // quotient ideal
857  tmp  = LR[2]; // varnames
858  s = size(LR[2]);
859  for (i=1; i<=D; i++)
860  {
861    for (j=1; j<=s; j++)
862    {
863      tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")";
864    }
865  }
866  for (i=1; i<=s; i++)
867  {
868    tmp[i] = string(tmp[i])+"("+string(1)+")";
869  }
870  L[2] = tmp;
871  list OrigNames = LR[2];
872  // ordering: d blocks of the ord on r
873  // try to get whether the ord on r is blockord itself
874  // TODO: make L(2) ordering! exponent is maximally 2
875  s = size(LR[3]);
876  if (s==2)
877  {
878    // not a blockord, 1 block + module ord
879    tmp = LR[3][s]; // module ord
880    for (i=1; i<=D; i++)
881    {
882      LR[3][s-1+i] = LR[3][1];
883    }
884    LR[3][s+D] = tmp;
885  }
886  if (s>2)
887  {
888    // there are s-1 blocks
889    int nb = s-1;
890    tmp = LR[3][s]; // module ord
891    for (i=1; i<=D; i++)
892    {
893      for (j=1; j<=nb; j++)
894      {
895        LR[3][i*nb+j] = LR[3][j];
896      }
897    }
898    //    size(LR[3]);
899    LR[3][nb*(D+1)+1] = tmp;
900  }
901  L[3] = LR[3];
902  def @R = ring(L);
903  //  setring @R;
904  return (@R);
905}
906example
907{
908  "EXAMPLE:"; echo = 2;
909  ring r = 0,(x,y,z),(dp(1),dp(2));
910  def A = freegbRing(2);
911  setring A;
912  A;
913}
914
915/* EXAMPLES:
916
917//static proc ex_shift()
918{
919  LIB "freegb.lib";
920  ring r = 0,(x,y,z),(dp(1),dp(2));
921  module M = [-1,x,y],[-7,y,y],[3,x,x];
922  module N = [1,x,y,x],[-1,y,x,y];
923  list L; L[1] = M; L[2] = N;
924  lst2str(L);
925  def U = crs(L,5);
926  setring U; U;
927  I;
928  poly p = I[2]; // I[8];
929  p;
930  system("stest",p,7,7,3); // error -> the world is ok
931  poly q1 = system("stest",p,1,7,3); //ok
932  poly q6 = system("stest",p,6,7,3); //ok
933  system("btest",p,3); //ok
934  system("btest",q1,3); //ok
935  system("btest",q6,3); //ok
936}
937
938//static proc test_shrink()
939{
940  LIB "freegb.lib";
941  ring r =0,(x,y,z),dp;
942  int d = 5;
943  def R = freegbRing(d);
944  setring R;
945  poly p1 = x(1)*y(2)*z(3);
946  poly p2 = x(1)*y(4)*z(5);
947  poly p3 = x(1)*y(1)*z(3);
948  poly p4 = x(1)*y(2)*z(2);
949  poly p5 = x(3)*z(5);
950  poly p6 = x(1)*y(1)*x(3)*z(5);
951  poly p7 = x(1)*y(2)*x(3)*y(4)*z(5);
952  poly p8 = p1+p2+p3+p4+p5 + p6 + p7;
953  p1; system("shrinktest",p1,3);
954  p2; system("shrinktest",p2,3);
955  p3; system("shrinktest",p3,3);
956  p4; system("shrinktest",p4,3);
957  p5; system("shrinktest",p5,3);
958  p6; system("shrinktest",p6,3);
959  p7; system("shrinktest",p7,3);
960  p8; system("shrinktest",p8,3);
961  poly p9 = p1 + 2*p2 + 5*p5 + 7*p7;
962  p9; system("shrinktest",p9,3);
963}
964
965//static proc ex2()
966{
967  option(prot);
968  LIB "freegb.lib";
969  ring r = 0,(x,y),dp;
970  module M = [-1,x,y],[3,x,x]; // 3x^2 - xy
971  def U = freegb(M,7);
972  lst2str(U);
973}
974
975//static proc ex_nonhomog()
976{
977  option(prot);
978  LIB "freegb.lib";
979  ring r = 0,(x,y,h),dp;
980  list L;
981  module M;
982  M = [-1,y,y],[1,x,x,x];  // x3-y2
983  L[1] = M;
984  M = [1,x,h],[-1,h,x];  // xh-hx
985  L[2] = M;
986  M = [1,y,h],[-1,h,y];  // yh-hy
987  L[3] = M;
988  def U = freegb(L,4);
989  lst2str(U);
990  // strange elements in the basis
991}
992
993//static proc ex_nonhomog_comm()
994{
995  option(prot);
996  LIB "freegb.lib";
997  ring r = 0,(x,y),dp;
998  module M = [-1,y,y],[1,x,x,x];
999  def U = freegb(M,5);
1000  lst2str(U);
1001}
1002
1003//static proc ex_nonhomog_h()
1004{
1005  option(prot);
1006  LIB "freegb.lib";
1007  ring r = 0,(x,y,h),(a(1,1),dp);
1008  module M = [-1,y,y,h],[1,x,x,x]; // x3 - y2h
1009  def U = freegb(M,6);
1010  lst2str(U);
1011}
1012
1013//static proc ex_nonhomog_h2()
1014{
1015  option(prot);
1016  LIB "freegb.lib";
1017  ring r = 0,(x,y,h),(dp);
1018  list L;
1019  module M;
1020  M = [-1,y,y,h],[1,x,x,x]; // x3 - y2h
1021  L[1] = M;
1022  M = [1,x,h],[-1,h,x]; // xh - hx
1023  L[2] = M;
1024  M = [1,y,h],[-1,h,y]; // yh - hy
1025  L[3] = M;
1026  def U = freegbasis(L,3);
1027  lst2str(U);
1029}
1030
1031
1032//static proc ex_nonhomog_3()
1033{
1034  option(prot);
1035  LIB "./freegb.lib";
1036  ring r = 0,(x,y,z),(dp);
1037  list L;
1038  module M;
1039  M = [1,z,y],[-1,x]; // zy - x
1040  L[1] = M;
1041  M = [1,z,x],[-1,y]; // zx - y
1042  L[2] = M;
1043  M = [1,y,x],[-1,z]; // yx - z
1044  L[3] = M;
1045  lst2str(L);
1046  list U = freegb(L,4);
1047  lst2str(U);
1049}
1050
1051//static proc ex_densep_2()
1052{
1053  option(prot);
1054  LIB "freegb.lib";
1055  ring r = (0,a,b,c),(x,y),(Dp); // deglex
1056  module M = [1,x,x], [a,x,y], [b,y,x], [c,y,y];
1057  lst2str(M);
1058  list U = freegb(M,5);
1059  lst2str(U);
1060  // a=b is important -> finite basis!!!
1061  module M = [1,x,x], [a,x,y], [a,y,x], [c,y,y];
1062  lst2str(M);
1063  list U = freegb(M,5);
1064  lst2str(U);
1065}
1066
1067// END COMMENTED EXAMPLES
1068
1069*/
1070
1071// 1. form a new ring
1072// 2. produce shifted generators
1073// 3. compute GB
1074// 4. skip shifted elts
1075// 5. go back to orig vars, produce strings/modules
1076// 6. return the result
1077
1078static proc freegbold(list LM, int d)
1079"USAGE:  freegbold(L, d);  L a list of modules, d an integer
1080RETURN:  ring
1081PURPOSE: compute the two-sided Groebner basis of an ideal, encoded by L in
1082the free associative algebra, up to degree d
1083EXAMPLE: example freegbold; shows examples
1084"
1085{
1086  // d = up to degree, will be shifted to d+1
1087  if (d<1) {"bad d"; return(0);}
1088
1089  int ppl = printlevel-voice+2;
1090  string err = "";
1091
1092  int i,j,s;
1093  def save = basering;
1094  // determine max no of places in the input
1095  int slm = size(LM); // numbers of polys in the ideal
1096  int sm;
1097  intvec iv;
1098  module M;
1099  for (i=1; i<=slm; i++)
1100  {
1101    // modules, e.g. free polynomials
1102    M  = LM[i];
1103    sm = ncols(M);
1104    for (j=1; j<=sm; j++)
1105    {
1106      //vectors, e.g. free monomials
1107      iv = iv, size(M[j])-1; // 1 place is reserved by the coeff
1108    }
1109  }
1110  int D = Max(iv); // max size of input words
1111  if (d<D) {"bad d"; return(LM);}
1112  D = D + d-1;
1113  //  D = d;
1114  list LR  = ringlist(save);
1115  list L, tmp;
1116  L[1] = LR[1]; // ground field
1117  L[4] = LR[4]; // quotient ideal
1118  tmp  = LR[2]; // varnames
1119  s = size(LR[2]);
1120  for (i=1; i<=D; i++)
1121  {
1122    for (j=1; j<=s; j++)
1123    {
1124      tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")";
1125    }
1126  }
1127  for (i=1; i<=s; i++)
1128  {
1129    tmp[i] = string(tmp[i])+"("+string(1)+")";
1130  }
1131  L[2] = tmp;
1132  list OrigNames = LR[2];
1133  // ordering: d blocks of the ord on r
1134  // try to get whether the ord on r is blockord itself
1135  // TODO: make L(2) ordering! exponent is maximally 2
1136  s = size(LR[3]);
1137  if (s==2)
1138  {
1139    // not a blockord, 1 block + module ord
1140    tmp = LR[3][s]; // module ord
1141    for (i=1; i<=D; i++)
1142    {
1143      LR[3][s-1+i] = LR[3][1];
1144    }
1145    LR[3][s+D] = tmp;
1146  }
1147  if (s>2)
1148  {
1149    // there are s-1 blocks
1150    int nb = s-1;
1151    tmp = LR[3][s]; // module ord
1152    for (i=1; i<=D; i++)
1153    {
1154      for (j=1; j<=nb; j++)
1155      {
1156        LR[3][i*nb+j] = LR[3][j];
1157      }
1158    }
1159    //    size(LR[3]);
1160    LR[3][nb*(D+1)+1] = tmp;
1161  }
1162  L[3] = LR[3];
1163  def @R = ring(L);
1164  setring @R;
1165  ideal I;
1166  poly @p;
1167  s = size(OrigNames);
1168  //  "s:";s;
1169  // convert LM to canonical vectors (no powers)
1170  setring save;
1171  kill M; // M was defined earlier
1172  module M;
1173  slm = size(LM); // numbers of polys in the ideal
1174  int sv,k,l;
1175  vector v;
1176  //  poly p;
1177  string sp;
1178  setring @R;
1179  poly @@p=0;
1180  setring save;
1181  for (l=1; l<=slm; l++)
1182  {
1183    // modules, e.g. free polynomials
1184    M  = LM[l];
1185    sm = ncols(M); // in intvec iv the sizes are stored
1186    for (i=0; i<=d-iv[l]; i++)
1187    {
1188      // modules, e.g. free polynomials
1189      for (j=1; j<=sm; j++)
1190      {
1191        //vectors, e.g. free monomials
1192        v  = M[j];
1193        sv = size(v);
1194        //        "sv:";sv;
1195        sp = "@@p = @@p + ";
1196        for (k=2; k<=sv; k++)
1197        {
1198          sp = sp + string(v[k])+"("+string(k-1+i)+")*";
1199        }
1200        sp = sp + string(v[1])+";"; // coef;
1201        setring @R;
1202        execute(sp);
1203        setring save;
1204      }
1205      setring @R;
1206      //      "@@p:"; @@p;
1207      I = I,@@p;
1208      @@p = 0;
1209      setring save;
1210    }
1211  }
1212  kill sp;
1213  // 3. compute GB
1214  setring @R;
1215  dbprint(ppl,"computing GB");
1216  //  ideal J = groebner(I);
1217  ideal J = slimgb(I);
1218  dbprint(ppl,J);
1219  // 4. skip shifted elts
1220  ideal K = select1(J,1..s); // s = size(OrigNames)
1221  dbprint(ppl,K);
1222  dbprint(ppl, "done with GB");
1223  // K contains vars x(1),...z(1) = images of originals
1224  // 5. go back to orig vars, produce strings/modules
1225  if (K[1] == 0)
1226  {
1227    "no reasonable output, GB gives 0";
1228    return(0);
1229  }
1230  int sk = size(K);
1231  int sp, sx, a, b;
1232  intvec x;
1233  poly p,q;
1234  poly pn;
1235  // vars in 'save'
1236  setring save;
1237  module N;
1238  list LN;
1239  vector V;
1240  poly pn;
1241  // test and skip exponents >=2
1242  setring @R;
1243  for(i=1; i<=sk; i++)
1244  {
1245    p  = K[i];
1246    while (p!=0)
1247    {
1249      //      "processing q:";q;
1251      sx = size(x);
1252      for(k=1; k<=sx; k++)
1253      {
1254        if ( x[k] >= 2 )
1255        {
1256          err = "skip: the value x[k] is " + string(x[k]);
1257          dbprint(ppl,err);
1258          //            return(0);
1259          K[i] = 0;
1260          p    = 0;
1261          q    = 0;
1262          break;
1263        }
1264      }
1265      p  = p - q;
1266    }
1267  }
1268  K  = simplify(K,2);
1269  sk = size(K);
1270  for(i=1; i<=sk; i++)
1271  {
1272    //    setring save;
1273    //    V  = 0;
1274    setring @R;
1275    p  = K[i];
1276    while (p!=0)
1277    {
1279      err =  "processing q:" + string(q);
1280      dbprint(ppl,err);
1282      sx = size(x);
1284      setring save;
1285      pn = imap(@R,pn);
1286      V  = V + leadcoef(pn)*gen(1);
1287      for(k=1; k<=sx; k++)
1288      {
1289        if (x[k] ==1)
1290        {
1291          a = k / s; // block number=a+1, a!=0
1292          b = k % s; // remainder
1293          //          printf("a: %s, b: %s",a,b);
1294          if (b == 0)
1295          {
1296            // that is it's the last var in the block
1297            b = s;
1298            a = a-1;
1299          }
1300          V = V + var(b)*gen(a+2);
1301        }
1302//         else
1303//         {
1304//           printf("error: the value x[k] is %s", x[k]);
1305//           return(0);
1306//         }
1307      }
1308      err = "V: " + string(V);
1309      dbprint(ppl,err);
1310      //      printf("V: %s", string(V));
1311      N = N,V;
1312      V  = 0;
1313      setring @R;
1314      p  = p - q;
1315      pn = 0;
1316    }
1317    setring save;
1318    LN[i] = simplify(N,2);
1319    N     = 0;
1320  }
1321  setring save;
1322  return(LN);
1323}
1324example
1325{
1326  "EXAMPLE:"; echo = 2;
1327  ring r = 0,(x,y,z),(dp(1),dp(2));
1328  module M = [-1,x,y],[-7,y,y],[3,x,x];
1329  module N = [1,x,y,x],[-1,y,x,y];
1330  list L; L[1] = M; L[2] = N;
1331  lst2str(L);
1332  def U = freegbold(L,5);
1333  lst2str(U);
1334}
1335
1336static proc sgb(ideal I, int d)
1337{
1338  // new code
1339  // map x_i to x_i(1) via map()
1340  //LIB "freegb.lib";
1341  def save = basering;
1342  //int d =7;// degree
1343  int nv = nvars(save);
1344  def R = freegbRing(d);
1345  setring R;
1346  int i;
1347  ideal Imap;
1348  for (i=1; i<=nv; i++)
1349  {
1350    Imap[i] = var(i);
1351  }
1352  //ideal I = x(1)*y(2), y(1)*x(2)+z(1)*z(2);
1353  ideal I = x(1)*x(2),x(1)*y(2) + z(1)*x(2);
1354  option(prot);
1355  //option(teach);
1356  ideal J = system("freegb",I,d,nv);
1357}
1358
1359static proc checkCeq()
1360{
1361  ring r = 0,(x,y),Dp;
1362  def A = freegbRing(4);
1363  setring A;
1364  A;
1365  // I = x2-xy
1366  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);
1367  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);
1368  ideal K = I,C;
1369  groebner(K);
1370}
1371
1372static proc exHom1()
1373{
1375  // z*y - x, z*x - y, y*x - z
1376  LIB "freegb.lib";
1377  LIB "elim.lib";
1378  ring r = 0,(x,y,z,h),dp;
1379  list L;
1380  module M;
1381  M = [1,z,y],[-1,x,h]; // zy - xh
1382  L[1] = M;
1383  M = [1,z,x],[-1,y,h]; // zx - yh
1384  L[2] = M;
1385  M = [1,y,x],[-1,z,h]; // yx - zh
1386  L[3] = M;
1387  lst2str(L);
1388  def U = crs(L,4);
1389  setring U;
1390  I = I,
1391    y(2)*h(3)+z(2)*x(3),     y(3)*h(4)+z(3)*x(4),
1392    y(2)*x(3)-z(2)*h(3),     y(3)*x(4)-z(3)*h(4);
1393  I = simplify(I,2);
1394  ring r2 = 0,(x(0..4),y(0..4),z(0..4),h(0..4)),dp;
1395  ideal J = imap(U,I);
1396  //  ideal K = homog(J,h);
1397  option(redSB);
1398  option(redTail);
1399  ideal L = groebner(J); //(K);
1400  ideal LL = sat(L,ideal(h))[1];
1401  ideal M = subst(LL,h,1);
1402  M = simplify(M,2);
1403  setring U;
1404  ideal M = imap(r2,M);
1405  lst2str(U);
1406}
1407
1408static proc test1()
1409{
1410  LIB "freegb.lib";
1411  ring r = 0,(x,y),Dp;
1412  int d = 10; // degree
1413  def R = freegbRing(d);
1414  setring R;
1415  ideal I = x(1)*x(2) - y(1)*y(2);
1416  option(prot);
1417  option(teach);
1418  ideal J = system("freegb",I,d,2);
1419  J;
1420}
1421
1422static proc test2()
1423{
1424  LIB "freegb.lib";
1425  ring r = 0,(x,y),Dp;
1426  int d = 10; // degree
1427  def R = freegbRing(d);
1428  setring R;
1429  ideal I = x(1)*x(2) - x(1)*y(2);
1430  option(prot);
1431  option(teach);
1432  ideal J = system("freegb",I,d,2);
1433  J;
1434}
1435
1436static proc test3()
1437{
1438  LIB "freegb.lib";
1439  ring r = 0,(x,y,z),dp;
1440  int d =5; // degree
1441  def R = freegbRing(d);
1442  setring R;
1443  ideal I = x(1)*y(2), y(1)*x(2)+z(1)*z(2);
1444  option(prot);
1445  option(teach);
1446  ideal J = system("freegb",I,d,3);
1447}
1448
1449static proc schur2-3()
1450{
1451  // nonhomog:
1452  //  h^4-10*h^2+9,f*e-e*f+h, h*2-e*h-2*e,h*f-f*h+2*f
1453  // homogenized with t
1454  //  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,
1455  // t*h - h*t, t*f - f*t, t*e - e*t
1456}
1457
1460RETURN:  ideal
1461ASSUME: there are at least i+j variables in the basering
1462PURPOSE: compute the ideal of Adem relations for i<2j in characteristic 0
1464"
1465{
1466  // produces Adem relations for i<2j in char 0
1467  // assume: 0<i<2j
1468  // requires presence of vars up to i+j
1469  if ( (i<0) || (i >= 2*j) )
1470  {
1471    ERROR("arguments out of range"); return(0);
1472  }
1473  ring @r = 0,(s(i+j..0)),lp;
1474  poly p,q;
1475  number n;
1476  int ii = i div 2; int k;
1477  // k=0 => s(0)=1
1478  n = binomial(j-1,i);
1479  q = n*s(i+j)*s(0);
1480  printf("k=0, term=%s",q);
1481  p = p + q;
1482  for (k=1; k<= ii; k++)
1483  {
1484    n = binomial(j-k-1,i-2*k);
1485    q = n*s(i+j-k)*s(k);;
1486    printf("k=%s, term=%s",k,q);
1487    p = p + q;
1488  }
1491  return(@r);
1492}
1493example
1494{
1495  "EXAMPLE:"; echo = 2;
1497  setring A;
1499}
1500
1501/*
15021,1: 0
15031,2: s(3)*s(0) == s(3) -> def for s(3):=s(1)s(2)
15052,2: s(3)*s(1) == s(1)s(2)s(1)
15061,3: 0 ( since 2*s(4)*s(0) = 0 mod 2)
15082,3: s(5)*s(0)+s(4)*s(1) == s(5)+s(4)*s(1)
15093,2: 0
15103,3: s(5)*s(1)
15111,4: 3*s(5)*s(0) == s(5)  -> def for s(5):=s(1)*s(4)
15132,4: 3*s(6)*s(0)+s(5)*s(1) == s(6) + s(5)*s(1) == s(6) + s(1)*s(4)*s(1)
15154,3: s(5)*s(2)
15163,4: s(7)*s(0)+2*s(6)*s(1) == s(7) -> def for s(7):=s(3)*s(4)
15174,4: s(7)*s(1)+s(6)*s(2)
1518*/
1519
1520/* s1,s2:
1521s1*s1 =0, s2*s2 = s1*s2*s1
1522*/
1523
1524/*
1525try char 0:
1526s1,s2:
1527s1*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)
1528hence 2==0! only in char 2
1529 */
1530
1531  // Adem rels modulo 2 are interesting
1532
1533//static
1534proc stringpoly2lplace(string s)
1535{
1536  // decomposes sentence into terms
1537  s = replace(s,newline,""); // get rid of newlines
1538  s = replace(s," ",""); // get rid of empties
1539  //arith symbols: +,-
1540  // decompose into words with coeffs
1541  list LS;
1542  int i,j,ie,je,k,cnt;
1543  // s[1]="-" situation
1544  if (s[1]=="-")
1545  {
1546    LS = stringpoly2lplace(string(s[2..size(s)]));
1547    LS[1] = string("-"+string(LS[1]));
1548    return(LS);
1549  }
1550  i = find(s,"-",2);
1551  // i==1 might happen if the 1st symbol coeff is negative
1552  j = find(s,"+");
1553  list LL;
1554  if (i==j)
1555  {
1556    "return a monomial";
1557    // that is both are 0 -> s is a monomial
1558    LS[1] = s;
1559    return(LS);
1560  }
1561  if (i==0)
1562  {
1563    "i==0 situation";
1564    // no minuses at all => pluses only
1565    cnt++;
1566    LS[cnt] = string(s[1..j-1]);
1567    s = s[j+1..size(s)];
1568    while (s!= "")
1569    {
1570      j = find(s,"+");
1571      cnt++;
1572      if (j==0)
1573      {
1574        LS[cnt] = string(s);
1575        s = "";
1576      }
1577      else
1578      {
1579        LS[cnt] = string(s[1..j-1]);
1580        s = s[j+1..size(s)];
1581      }
1582    }
1583    return(LS);
1584  }
1585  if (j==0)
1586  {
1587    "j==0 situation";
1588    // no pluses at all except the lead coef => the rest are minuses only
1589    cnt++;
1590    LS[cnt] = string(s[1..i-1]);
1591    s = s[i..size(s)];
1592    while (s!= "")
1593    {
1594      i = find(s,"-",2);
1595      cnt++;
1596      if (i==0)
1597      {
1598        LS[cnt] = string(s);
1599        s = "";
1600      }
1601      else
1602      {
1603        LS[cnt] = string(s[1..i-1]);
1604        s = s[i..size(s)];
1605      }
1606    }
1607    return(LS);
1608  }
1609  // now i, j are nonzero
1610  if (i>j)
1611  {
1612    "i>j situation";
1613    // + comes first, at place j
1614    cnt++;
1615    //    "cnt:"; cnt; "j:"; j;
1616    LS[cnt] = string(s[1..j-1]);
1617    s = s[j+1..size(s)];
1618    LL = stringpoly2lplace(s);
1619    LS = LS + LL;
1620    kill LL;
1621    return(LS);
1622  }
1623  else
1624  {
1625    "j>i situation";
1626    // - might come first, at place i
1627    if (i>1)
1628    {
1629      cnt++;
1630      LS[cnt] = string(s[1..i-1]);
1631      s = s[i..size(s)];
1632    }
1633    else
1634    {
1635      // i==1->  minus at leadcoef
1636      ie = find(s,"-",i+1);
1637      je = find(s,"+",i+1);
1638      if (je == ie)
1639      {
1640         "ie=je situation";
1641        //monomial
1642        cnt++;
1643        LS[cnt] = s;
1644        return(LS);
1645      }
1646      if (je < ie)
1647      {
1648         "je<ie situation";
1649        // + comes first
1650        cnt++;
1651        LS[cnt] = s[1..je-1];
1652        s = s[je+1..size(s)];
1653      }
1654      else
1655      {
1656        // ie < je
1657         "ie<je situation";
1658        cnt++;
1659        LS[cnt] = s[1..ie-1];
1660        s = s[ie..size(s)];
1661      }
1662    }
1663    "going into recursion with "+s;
1664    LL = stringpoly2lplace(s);
1665    LS = LS + LL;
1666    return(LS);
1667  }
1668}
1669example
1670{
1671  "EXAMPLE:"; echo = 2;
1672  string s = "x*y+y*z+z*t"; // + only
1673  stringpoly2lplace(s);
1674  string s2 = "x*y - y*z-z*t*w*w"; // +1, - only
1675  stringpoly2lplace(s2);
1676  string s3 = "-x*y + y - 2*x +7*w*w*w";
1677  stringpoly2lplace(s3);
1678}
1679
1681{
1682  // adds places to the list of strings
1683  // according to their order in the list
1684  int s = size(L);
1685  int i;
1686  for (i=1; i<=s; i++)
1687  {
1688    if (typeof(L[i]) == "string")
1689    {
1690      L[i] = L[i] + "(" + string(i) + ")";
1691    }
1692    else
1693    {
1694      ERROR("entry of type string expected");
1695      return(0);
1696    }
1697  }
1698  return(L);
1699}
1700example
1701{
1702  "EXAMPLE:"; echo = 2;
1703  string a = "f1";   string b = "f2";
1704  list L = a,b,a;
1706}
1707
1708//static
1709proc sent2lplace(string s)
1710{
1711  // SENTence of words TO LetterPLACE
1712  list L =   stringpoly2lplace(s);
1713  int i; int ss = size(L);
1714  for(i=1; i<=ss; i++)
1715  {
1716    L[i] = str2lplace(L[i]);
1717  }
1718  return(L);
1719}
1720example
1721{
1722  "EXAMPLE:"; echo = 2;
1723  ring r = 0,(f2,f1),dp;
1724  string s = "f2*f1*f1 - 2*f1*f2*f1+ f1*f1*f2";
1725  sent2lplace(s);
1726}
1727
1728static proc testnumber(string s)
1729{
1730  string t;
1731  if (s[1]=="-")
1732  {
1733    // two situations: either there's a negative number
1734    t = s[2..size(s)];
1735    if (testnumber(t))
1736    {
1737      //a negative number
1738    }
1739    else
1740    {
1741      // a variable times -1
1742    }
1743    // or just a "-" for -1
1744  }
1745  t = "ring @r=(";
1746  t = t + charstr(basering)+"),";
1747  t = t + string(var(1))+",dp;";
1748  //  write(":w tstnum.tst",t);
1749  t = t+ "number @@Nn = " + s + ";"+"\$";
1750  write(":w tstnum.tst",t);
1751  string runsing = system("Singular");
1752  int k;
1753  t = runsing+ " -teq <tstnum.tst >tstnum.out";
1754  k = system("sh",t);
1755  if (k!=0)
1756  {
1757    ERROR("Problems running Singular");
1758  }
1759  int i = system("sh", "grep error tstnum.out > /dev/NULL");
1760  if (i!=0)
1761  {
1762    // no error: s is a number
1763    i = 1;
1764  }
1765  k = system("sh","rm tstnum.tst tstnum.out > /dev/NULL");
1766  return(i);
1767}
1768example
1769{
1770  "EXAMPLE:"; echo = 2;
1771  ring r = (0,a),x,dp;
1772  string s = "a^2+7*a-2";
1773  testnumber(s);
1774  s = "b+a";
1775  testnumber(s);
1776}
1777
1778proc str2lplace(string s)
1779{
1780  // converts a word (monomial) with coeff into letter-place
1781  // string: coef*var1^exp1*var2^exp2*...varN^expN
1782  s = strpower2rep(s); // expand powers
1783  if (size(s)==0) { return(0); }
1784  int i,j,k,insC;
1785  string a,b,c,d,t;
1786  // 1. get coeff
1787  i = find(s,"*");
1788  if (i==0) { return(s); }
1789  list VN;
1790  c = s[1..i-1]; // incl. the case like (-a^2+1)
1791  int tn = testnumber(c);
1792  if (tn == 0)
1793  {
1794    // failed test
1795    if (c[1]=="-")
1796    {
1797      // two situations: either there's a negative number
1798      t = c[2..size(c)];
1799      if (testnumber(t))
1800      {
1801         //a negative number
1802        // nop here
1803      }
1804      else
1805      {
1806         // a variable times -1
1807          c = "-1";
1808          j++; VN[j] = t; //string(c[2..size(c)]);
1809          insC = 1;
1810      }
1811    }
1812    else
1813    {
1814      // just a variable with coeff 1
1815          j++; VN[j] = string(c);
1816          c = "1";
1817          insC = 1;
1818    }
1819  }
1820 // get vars
1821  t = s;
1822  //  t = s[i+1..size(s)];
1823  k = size(t); //j = 0;
1824  while (k>0)
1825  {
1826    t = t[i+1..size(t)]; //next part
1827    i = find(t,"*"); // next *
1828    if (i==0)
1829    {
1830      // last monomial
1831      j++;
1832      VN[j] = t;
1833      k = size(t);
1834      break;
1835    }
1836    b = t[1..i-1];
1837    //    print(b);
1838    j++;
1839    VN[j] = b;
1840    k = size(t);
1841  }
1843  VN[size(VN)+1] = string(c);
1844  return(VN);
1845}
1846example
1847{
1848  "EXAMPLE:"; echo = 2;
1849  ring r = (0,a),(f2,f1),dp;
1850  str2lplace("-2*f2^2*f1^2*f2");
1851  str2lplace("-f1*f2");
1852  str2lplace("(-a^2+7a)*f1*f2");
1853}
1854
1855static proc strpower2rep(string s)
1856{
1857  // makes x*x*x*x out of x^4 ., rep statys for repetitions
1858  // looks for "-" problem
1859  // exception: "-" as coeff
1860  string ex,t;
1861  int i,j,k;
1862
1863  i = find(s,"^"); // first ^
1864  if (i==0) { return(s); } // no ^ signs
1865
1866  if (s[1] == "-")
1867  {
1868    // either -coef or -1
1869    // got the coeff:
1870    i = find(s,"*");
1871    if (i==0)
1872    {
1873      // no *'s   => coef == -1 or s == -23
1874      i = size(s)+1;
1875    }
1876    t = string(s[2..i-1]); // without "-"
1877    if ( testnumber(t) )
1878    {
1879      // a good number
1880      t = strpower2rep(string(s[2..size(s)]));
1881      t = "-"+t;
1882      return(t);
1883    }
1884    else
1885    {
1886      // a variable
1887      t = strpower2rep(string(s[2..size(s)]));
1888      t = "-1*"+ t;
1889      return(t);
1890    }
1891  }
1892  // the case when leadcoef is a number in ()
1893  if (s[1] == "(")
1894  {
1895    i = find(s,")",2);    // must be nonzero
1896    t = s[2..i-1];
1897    if ( testnumber(t) )
1898    {
1899      // a good number
1900    }
1901    else {"strpower2rep: bad number as coef";}
1902    ex = string(s[i+2..size(s)]); // 2 because of *
1903    ex =  strpower2rep(ex);
1904    t = "("+t+")*"+ex;
1905    return(t);
1906  }
1907
1908  i = find(s,"^"); // first ^
1909  j = find(s,"*",i+1); // next * == end of ^
1910  if (j==0)
1911  {
1912    ex = s[i+1..size(s)];
1913  }
1914  else
1915  {
1916    ex = s[i+1..j-1];
1917  }
1918  execute("int @exp = " + ex + ";"); //@exp = exponent
1919  // got varname
1920  for (k=i-1; k>0; k--)
1921  {
1922    if (s[k] == "*") break;
1923  }
1924  string varn = s[k+1..i-1];
1925  //  "varn:";  varn;
1926  string pref;
1927  if (k>0)
1928  {
1929    pref = s[1..k]; // with * on the k-th place
1930  }
1931  //  "pref:";  pref;
1932  string suf;
1933  if ( (j>0) && (j+1 <= size(s)) )
1934  {
1935    suf = s[j+1..size(s)]; // without * on the 1st place
1936  }
1937  //  "suf:"; suf;
1938  string toins;
1939  for (k=1; k<=@exp; k++)
1940  {
1941    toins = toins + varn+"*";
1942  }
1943  //  "toins: ";  toins;
1944  if (size(suf) == 0)
1945  {
1946    toins = toins[1..size(toins)-1]; // get rid of trailing *
1947  }
1948  else
1949  {
1950    suf = strpower2rep(suf);
1951  }
1952  ex = pref + toins + suf;
1953  return(ex);
1954  //  return(strpower2rep(ex));
1955}
1956example
1957{
1958  "EXAMPLE:"; echo = 2;
1959  ring r = (0,a),(x,y,z,t),dp;
1960  strpower2rep("-x^4");
1961  strpower2rep("-2*x^4*y^3*z*t^2");
1962  strpower2rep("-a^2*x^4");
1963}
1964
1965proc Liebr(poly a, poly b, list #)
1966"USAGE:  Liebr(a,b[,N]); a,b letterplace polynomials, N an optional integer
1967RETURN:  poly
1968PURPOSE: compute the Lie bracket [a,b] = ab - ba between letterplace polynomials
1969NOTE: if N>1 is specified, then the left normed bracket [a,[...[a,b]]]] is computed.
1970EXAMPLE: example Liebr; shows examples
1971"
1972{
1973  // alias ppLiebr;
1974  //if int N is given compute [a,[...[a,b]]]] left normed bracket
1975  poly q;
1976  int N=1;
1977  if (size(#)>0)
1978  {
1979    if (typeof(#[1])=="int")
1980    {
1981      N = int(#[1]);
1982    }
1983  }
1984  if (N<=0) { return(q); }
1985  while (b!=0)
1986  {
1987    q = q + pmLiebr(a,lead(b));
1988    b = b - lead(b);
1989  }
1990  int i;
1991  if (N >1)
1992  {
1993    for(i=1; i<=N; i++)
1994    {
1995      q = Liebr(a,q);
1996    }
1997  }
1998  return(q);
1999}
2000example
2001{
2002  "EXAMPLE:"; echo = 2;
2003  ring r = 0,(x(1),y(1),x(2),y(2),x(3),y(3),x(4),y(4)),dp;
2004  poly a = x(1)*y(2); poly b = y(1);
2005  int uptodeg=4; int lV=2;
2006  export uptodeg; export lV;
2007  Liebr(a,b);
2008  Liebr(x(1),y(1),2);
2009}
2010
2011static proc pmLiebr(poly a, poly b)
2012{
2013  //  a poly, b mono
2014  poly s;
2015  while (a!=0)
2016  {
2018    a = a - lead(a);
2019  }
2020  return(s);
2021}
2022
2023//proc pshift(poly a, int i, int uptodeg, int lV)
2024static proc pshift(poly a, int i)
2025{
2026  // shifts a monomial a by i
2027  // calls pLPshift(p,sh,uptodeg,lVblock);
2028  return(system("stest",a,i,uptodeg,lV));
2029}
2030
2031static proc mmLiebr(poly a, poly b)
2032{
2033  // a,b, monomials
2036  int sa = deg(a);
2037  int sb = deg(b);
2038  poly v = a*pshift(b,sa) - b*pshift(a,sb);
2039  return(v);
2040}
2041
2042static proc test_shift()
2043{
2044  LIB "freegb.lib";
2045  ring r = 0,(a,b),dp;
2046  int d =5;
2047  def R = freegbRing(d);
2048  setring R;
2049  int uptodeg = d; export uptodeg;
2050  int lV = 2; export lV;
2051  poly p = mmLiebr(a(1),b(1));
2052  poly p = Liebr(a(1),b(1));
2053}
2054
2055proc Serre(intmat A, int zu)
2056"USAGE:  Serre(A,z); A an intmat, z an int
2057RETURN:  ideal
2058ASSUME: basering has a letterplace ring structure and
2059@*    A is a generalized Cartan matrix with integer entries
2060PURPOSE: compute the ideal of Serre's relations associated to A
2061EXAMPLE: example Serre; shows examples
2062"
2063{
2064  // zu = 1 -> with commutators [f_i,f_j]; zu == 0 without them
2065  // suppose that A is cartan matrix
2066  // then Serre's relations are
2067  // (ad f_j)^{1-A_{ij}} ( f_i)
2068  int ppl = printlevel-voice+2;
2069  int n = ncols(A); // hence n variables
2070  int i,j,k,el;
2071  poly p,q;
2072  ideal I;
2073  for (i=1; i<=n; i++)
2074  {
2075    for (j=1; j<=n; j++)
2076    {
2077      el = 1 - A[i,j];
2078      //     printf("i:%s, j: %s, l: %s",i,j,l);
2079      dbprint(ppl,"i, j, l: ",i,j,el);
2080      //      if ((i!=j) && (l >0))
2081      //      if ( (i!=j) &&  ( ((zu ==0) &&  (l >=2)) || ((zu ==1) &&  (l >=1)) ) )
2082      if ((i!=j) && (el >0))
2083      {
2084        q = Liebr(var(j),var(i));
2085        dbprint(ppl,"first bracket: ",q);
2086        //        if (l >=2)
2087        //        {
2088          for (k=1; k<=el-1; k++)
2089          {
2090            q = Liebr(var(j),q);
2091            dbprint(ppl,"further bracket:",q);
2092          }
2093          //        }
2094      }
2095      if (q!=0) { I = I,q; q=0;}
2096    }
2097  }
2098  I = simplify(I,2);
2099  return(I);
2100}
2101example
2102{
2103  "EXAMPLE:"; echo = 2;
2104  intmat A[3][3] =
2105    2, -1, 0,
2106    -1, 2, -3,
2107    0, -1, 2; // G^1_2 Cartan matrix
2108  ring r = 0,(f1,f2,f3),dp;
2109  int uptodeg = 5; int lV = 3;
2110  export uptodeg; export lV;
2111  def R = freegbRing(uptodeg);
2112  setring R;
2113  ideal I = Serre(A,1); I = simplify(I,1+2+8);
2114  I;
2115}
2116
2117/* setup for older example:
2118  intmat A[2][2] = 2, -1, -1, 2; // sl_3 == A_2
2119  ring r = 0,(f1,f2),dp;
2120  int uptodeg = 5; int lV = 2;
2121*/
2122
2123proc lp2lstr(ideal K, def save)
2124"USAGE:  lp2lstr(K,s); K an ideal, s a ring
2125RETURN:  nothing (exports object LN into s)
2126ASSUME: basering has a letterplace ring structure
2127PURPOSE: converts letterplace ideal to list of modules
2128NOTE: useful as preprocessing to 'lst2str'
2129EXAMPLE: example lp2lstr; shows examples
2130"
2131{
2132  def @R = basering;
2133  string err;
2134  int s = nvars(save);
2135  int i,j,k;
2136    // K contains vars x(1),...z(1) = images of originals
2137  // 5. go back to orig vars, produce strings/modules
2138  int sk = size(K);
2139  int sp, sx, a, b;
2140  intvec x;
2141  poly p,q;
2142  poly pn;
2143  // vars in 'save'
2144  setring save;
2145  module N;
2146  list LN;
2147  vector V;
2148  poly pn;
2149  // test and skip exponents >=2
2150  setring @R;
2151  for(i=1; i<=sk; i++)
2152  {
2153    p  = K[i];
2154    while (p!=0)
2155    {
2157      //      "processing q:";q;
2159      sx = size(x);
2160      for(k=1; k<=sx; k++)
2161      {
2162        if ( x[k] >= 2 )
2163        {
2164          err = "skip: the value x[k] is " + string(x[k]);
2165          dbprint(ppl,err);
2166          //            return(0);
2167          K[i] = 0;
2168          p    = 0;
2169          q    = 0;
2170          break;
2171        }
2172      }
2173      p  = p - q;
2174    }
2175  }
2176  K  = simplify(K,2);
2177  sk = size(K);
2178  for(i=1; i<=sk; i++)
2179  {
2180    //    setring save;
2181    //    V  = 0;
2182    setring @R;
2183    p  = K[i];
2184    while (p!=0)
2185    {
2187      err =  "processing q:" + string(q);
2188      dbprint(ppl,err);
2190      sx = size(x);
2192      setring save;
2193      pn = imap(@R,pn);
2194      V  = V + leadcoef(pn)*gen(1);
2195      for(k=1; k<=sx; k++)
2196      {
2197        if (x[k] ==1)
2198        {
2199          a = k / s; // block number=a+1, a!=0
2200          b = k % s; // remainder
2201          //          printf("a: %s, b: %s",a,b);
2202          if (b == 0)
2203          {
2204            // that is it's the last var in the block
2205            b = s;
2206            a = a-1;
2207          }
2208          V = V + var(b)*gen(a+2);
2209        }
2210      }
2211      err = "V: " + string(V);
2212      dbprint(ppl,err);
2213      //      printf("V: %s", string(V));
2214      N = N,V;
2215      V  = 0;
2216      setring @R;
2217      p  = p - q;
2218      pn = 0;
2219    }
2220    setring save;
2221    LN[i] = simplify(N,2);
2222    N     = 0;
2223  }
2224  setring save;
2225  export LN;
2226  //  return(LN);
2227}
2228example
2229{
2230  "EXAMPLE:"; echo = 2;
2231  intmat A[2][2] = 2, -1, -1, 2; // sl_3 == A_2
2232  ring r = 0,(f1,f2),dp;
2233  int uptodeg = 3; int lV = 2;
2234  export uptodeg; export lV;
2235  def R = freegbRing(uptodeg);
2236  setring R;
2237  ideal I = Serre(A,1);
2238  lp2lstr(I,r);
2239  setring r;
2240  lst2str(LN,1);
2241  kill uptodeg; kill lV;
2242}
2243
2244proc strList2poly(list L)
2245{
2246  //  list L comes from sent2lplace (which takes a poly on the input)
2247  // each entry of L is a sublist with the coef on the last place
2248  int s = size(L); int t;
2249  int i,j;
2250  list M;
2251  poly p,q;
2252  string Q;
2253  for(i=1; i<=s; i++)
2254  {
2255    M = L[i];
2256    t = size(M);
2257    //    q = M[t]; // a constant
2258    Q = string(M[t]);
2259    for(j=1; j<t; j++)
2260    {
2261      //      q = q*M[j];
2262      Q = Q+"*"+string(M[j]);
2263    }
2264    execute("q="+Q+";");
2265    //    q;
2266    p = p + q;
2267  }
2268  kill Q;
2269  return(p);
2270}
2271example
2272{
2273  "EXAMPLE:"; echo = 2;
2274  ring r =0,(x,y,z,t),Dp;
2275  def A = freegbRing(4);
2276  setring A;
2277  string t = "-2*y*z*y*z + y*t*z*z - z*x*x*y  + 2*z*y*z*y";
2278  list L = sent2lplace(t);
2279  L;
2280  poly p = strList2poly(L);
2281  p;
2282}
2283
2284proc file2lplace(string fname)
2285"USAGE:  file2lplace(fnm);  fnm a string
2286RETURN:  ideal
2287PURPOSE: convert the contents of the file fnm into ideal of polynomials in free algebra
2288EXAMPLE: example file2lplace; shows examples
2289"
2290{
2291  // format: from the usual string to letterplace
2293  // assume: file is a comma-sep list of polys
2294  // the vars are declared before
2295  // the file ends with ";"
2296  string t; int i;
2297  ideal I;
2298  list tst;
2299  while (s!="")
2300  {
2301    i = find(s,",");
2302    "i"; i;
2303    if (i==0)
2304    {
2305      i = find(s,";");
2306      if (i==0)
2307      {
2308        // no ; ??
2309         "no colon or semicolon found anymore";
2310         return(I);
2311      }
2312      // no "," but ";" on the i-th place
2313      t = s[1..i-1];
2314      s = "";
2315      "processing: "; t;
2316      tst = sent2lplace(t);
2317      tst;
2318      I = I, strList2poly(tst);
2319      return(I);
2320    }
2321    // here i !=0
2322    t = s[1..i-1];
2323    s = s[i+1..size(s)];
2324    "processing: "; t;
2325    tst = sent2lplace(t);
2326    tst;
2327    I = I, strList2poly(tst);
2328  }
2329  return(I);
2330}
2331example
2332{
2333  "EXAMPLE:"; echo = 2;
2334  ring r =0,(x,y,z,t),dp;
2335  def A = freegbRing(4);
2336  setring A;
2337  string fn = "myfile";
2338  string s1 = "z*y*y*y - 3*y*z*x*y  + 3*y*y*z*y - y*x*y*z,";
2339  string s2 = "-2*y*x*y*z + y*y*z*z - z*z*y*y + 2*z*y*z*y,";
2340  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;";
2341  write(":w "+fn,s1);  write(":a "+fn,s2);   write(":a "+fn,s3);
2343  ideal I = file2lplace(fn);
2344  I;
2345}
2346
2347/* EXAMPLES AGAIN:
2348//static proc get_ls3nilp()
2349{
2350//first app of file2lplace
2351  ring r =0,(x,y,z,t),dp;
2352  int d = 10;
2353  def A = freegbRing(d);
2354  setring A;
2355  ideal I = file2lplace("./ls3nilp.bg");
2356  // and now test the correctness: go back from lplace to strings
2357  lp2lstr(I,r);
2358  setring r;
2359  lst2str(LN,1); // agree!
2360}
2361
2362//static proc doc_example()
2363{
2364  LIB "freegb.lib";
2365  ring r = 0,(x,y,z),dp;
2366  int d =4; // degree bound
2367  def R = freegbRing(d);
2368  setring R;
2369  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);
2370  option(redSB);option(redTail);
2371  ideal J = system("freegb",I,d,nvars(r));
2372  J;
2373  // visualization:
2374  lp2lstr(J,r); // export an object called LN to the ring r
2375  setring r;  // change to the ring r
2376  lst2str(LN,1); // output the strings
2377}
2378
2379*/
2380
2381// TODO:
2382// multiply two letterplace polynomials, lpMult
2383// reduction/ Normalform? needs kernel stuff
2384
Note: See TracBrowser for help on using the repository browser.