# Ticket #54: freegb.lib.merge

File freegb.lib.merge, 55.0 KB (added by Oleksandr , 14 years ago)

Merge of new library together with my comments/corrections

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