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

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