# source:git/Singular/LIB/sagbi.lib@043cba

spielwiese
Last change on this file since 043cba was 043cba, checked in by Hans Schönemann <hannes@…>, 16 years ago
• Property mode set to `100644`
File size: 16.2 KB
Line
1//////////////////////////////////////////////////////////////////////////////
2version="\$Id: sagbi.lib,v 1.11 2007-07-13 12:31:06 Singular Exp \$";
3category="Commutative Algebra";
4info="
5LIBRARY:  sagbi.lib  Compute Subalgebras bases Analogous to Groebner bases for ideals
6AUTHORS: Gerhard Pfister,        pfister@mathematik.uni-kl.de,
7@*       Anen Lakhal,            alakhal@mathematik.uni-kl.de
8
9PROCEDURES:
10 sagbiRreduction(p,I);  perform one step subalgebra reducton (for short S-reduction) of p w.r.t I
11 sagbiSPoly(I);   compute the S-polynomials of the Subalgebra defined by the genartors of I
12 sagbiNF(id,I);   perform iterated S-reductions in order to compute Subalgebras normal forms
13 sagbi(I);        construct SAGBI basis for the Subalgebra defined by I
14 sagbiPart(I);    construct partial SAGBI basis for the Subalgebra defined by I
15";
16
17LIB "algebra.lib";
18LIB "elim.lib";
19
20///////////////////////////////////////////////////////////////////////////////
21proc sagbiSPoly(id ,list #)
22"USAGE: sagbiSPoly(id [,n]); id ideal, n positive integer.
23RETURN: an ideal
24@format
25      - If (n=0 or default) an ideal, whose generators are the S-polynomials.
26      - If (n=1) a  list of size 2:
27        the first element of this list is the ideal of S-polynomials.
28        the second element of this list is the ring in which is defined
29        the ideal of algebraic relations.
30@end format
31EXAMPLE: example sagbiSPoly; show an example "
32{
33  if(size(#)==0)
34  {
35    #[1]=0;
36  }
37  degBound=0;
38  def bsr=basering;
39  ideal vars=maxideal(1);
40  ideal B=ideal(bsr);//when the basering is quotient ring this "type casting"
41                    //gives th quotient ideal.
42  int b=size(B);
43
44  //In quotient rings,SINGULAR does not reduce polynomials w.r.t the
45  //quotient ideal,therefore we should first 'reduce';if it is necessary for
46  //computations to have a uniquely determined representant for each equivalent
47  //class,which is the case of this procedure.
48
49  if(b!=0)
50  {
51    id =reduce(id,groebner(0));
52  }
53  int n,m=nvars(bsr),ncols(id);
54  int z;
55  string mp=string(minpoly);
56  ideal P;
57  list L;
58
59  if(id==0)
60  {
61    if(#[1]==0)
62    {
63      return(P);
64    }
65    else
66    {
67      return(L);
68    }
69  }
70  else
71  {
72    //==================create anew ring with extra variables================
73
74    execute("ring R1=("+charstr(bsr)+"),("+varstr(bsr)+",@y(1..m)),(dp(n),dp(m));");
75    execute("minpoly=number("+mp+");");
76    ideal id=imap(bsr,id);
77    ideal A;
78
79    for(z=1;z<=m;z++)
80    {
82    }
83
84    A=groebner(A);
85    ideal kern=nselect(A,1,n);// "kern" is the kernel of the ring map:
87                        //"kern" is the ideal of algebraic relations between
89    export kern,A;// we export:
90                  // * the ideal A  to avoid useless computations
91                  //   between 2 steps in sagbi procedure.
92                  // * the ideal kern : some times we can get intresting
93                  //   informations from this ideal.
94    setring bsr;
95    map phi=R1,vars,id;
96
97    // the sagbiSPolynomials are the image by phi of the generators of kern
98
99    P=simplify(phi(kern),1);
100    if(#[1]==0)
101    {
102      return(P);
103    }
104    else
105    {
106      L=P,R1;
107      kill phi,vars;
108
109      dbprint(printlevel-voice+3,"
110// 'sagbiSPoly' created a ring as 2nd element of the list.
111// The ring contains the ideal 'kern'  of algebraic relations between the
112//leading terms of the generators of I.
113// To access to this ring and see 'kern' you should give the ring a name,
114// e.g.:
115               def S = L[2]; setring S; kern;
116      ");
117      return(L);
118    }
119  }
120}
121example
122{ "EXAMPLE:"; echo = 2;
123 ring r=0, (x,y),dp;
124 poly f1,f2,f3,f4=x2,y2,xy+y,2xy2;
125 ideal I=f1,f2,f3,f4;
126 sagbiSPoly(I);
127 list L=sagbiSPoly(I,1);
128 L[1];
129 def S= L[2]; setring S; kern;
130}
131
132///////////////////////////////////////////////////////////////////////////////
133static proc std1(ideal J,ideal I,list #)
134  // I is contained in J, and it is assumed to be a standard bases!
135  //This procedure computes a Standard basis for J from I one's
136  //This procedure is essential for Spoly1 procedure.
137{
138  def br=basering;
139  int tt;
140  ideal Res,@result;
141
142  if(size(#)>0) {tt=#[1];}
143
144  if(size(I)==0) {@result=groebner(J);}
145
146  if((size(I)!=0) && (size(J)-size(I)>=1))
147  {
148    qring Q=I;
149    ideal J=fetch(br,J);
150    J=groebner(J);
151    setring br;
152    Res=fetch(Q,J);// Res contains the generators that we add to I
153                   // to get the generators of std(J);
154    @result=Res+I;
155  }
156
157  if(tt==0) { return(@result);}
158  else      { return(Res);}
159}
160
161///////////////////////////////////////////////////////////////////////////////
162
163static proc Spoly1(list l,ideal I,ideal J,int a)
164  //an implementation of SAGBI construction Algorithm using Spoly
165  //procedure leads to useless computations and affect the efficiency
166  //of SAGBI bases computations. This procedure is a variant of Spoly
167  //in order to avoid these useless compuations.
168{
169  degBound=0;
170  def br=basering;
171  ideal vars=maxideal(1);
172  ideal B=ideal(br);
173  int b=size(B);
174
175  if(b!=0)
176  {
177    I=reduce(I,groebner(0));
178    J=reduce(J,groebner(0));
179  }
180  int n,ii,jj=nvars(br),ncols(I),ncols(J);
181  int z;
182  list @L;
183  string mp =string(minpoly);
184
185  if(size(J)==0)
186  {
187    @L =sagbiSPoly(I,1);
188  }
189  else
190  {
191    ideal @sum=I+J;
192    ideal P1;
193    ideal P=l[1];//P is the ideal of spolynomials of I;
194    def R=l[2];setring R;int kk=nvars(R);
195    ideal J=fetch(br,J);
196
197    //================create a new ring with extra variables==============
198    execute("ring R1=("+charstr(R)+"),("+varstr(R)+",@y((ii+1)..(ii+jj))),(dp(n),dp(kk+jj-n));");
199    // *levandov: would it not be easier and better to use
200    // ring @Y = char(R),(@y((ii+1)..(ii+jj))),dp;
201    // def R1 = R + @Y;
202    // setring R1;
203    // -> thus
204    ideal kern1;
205    ideal A=fetch(R,A);
206    attrib(A,"isSB",1);
207    ideal J=fetch(R,J);
208    ideal kern=fetch(R,kern);
209    ideal A1;
210    for(z=1;z<=jj;z++)
211    {
213    }
214    A1=A+A1;
215    ideal @Res=std1(A1,A,1);// the generators of @Res are whose we have to add
216                            // to A to get std(A1).
217    A=A+@Res;
218    kern1=nselect(@Res,1,n);
219    kern=kern+kern1;
220    export kern,kern1,A;
221    setring br;
222    map phi=R1,vars,@sum;
223    P1=simplify(phi(kern1),1);//P1 is th ideal we add to P to get the ideal
224                              //of Spolynomials of @sum.
225    P=P+P1;
226
227    if (a==1)
228    {
229      @L=P,R1;
230      kill phi,vars;
231      dbprint(printlevel-voice+3,"
232// 'Spoly1' created a ring as 2nd element of the list.
233// The ring contains the ideal 'kern'  of algebraic relations between the
234//generators of I+J.
235// To access to this ring and see 'kern' you should give the ring a name,
236// e.g.:
237               def @ring = L[2]; setring @ring ; kern;
238      ");
239    }
240    if(a==2)
241    {
242      @L=P1,R1;
243      kill phi,vars;
244    }
245  }
246  return(@L);
247}
248///////////////////////////////////////////////////////////////////////////////
249
250proc sagbiReduction(poly p,ideal dom,list #)
251"USAGE: sagbiReduction(p,dom[,n]); p poly , dom  ideal
252RETURN: a polynomial, after one step subalgebra reduction
253@format
254    Three algorithm variants are used to perform subalgebra reduction.
255    The positive interger n determines which variant should be used.
256    n may take the values 0 (default), 1 or 2.
257@end format
258EXAMPLE: sagbiReduction; shows an example"
259{
260  def bsr=basering;
261  ideal B=ideal(bsr);//When the basering is quotient ring  this type casting
262                     // gives the quotient ideal.
263  int b=size(B);
264  int n=nvars(bsr);
265
266  //In quotient rings, SINGULAR, usually does not reduce polynomials w.r.t the
267  //quotient ideal,therefore we should first  reduce ,when it is necessary for computations,
268  // to have a uniquely determined representant for each equivalent
269  //class,which is the case of this algorithm.
270
271  if(b !=0) //means that the basering is a quotient ring
272  {
273    p=reduce(p,std(0));
274    dom=reduce(dom,std(0));
275  }
276
277  int i,choose;
278  int z=ncols(dom);
279
280  if((size(#)>0) && (typeof(#[1])=="int"))
281  {
282    choose = #[1];
283  }
284  if (size(#)>1)
285  {
286    choose =#[2];
287  }
288
289  //=======================first algorithm(default)=========================
290  if ( choose == 0 )
291  {
293    if( L[1]==1 )
294    {
295      // the ring L[2] = char(bsr),(x(1..nvars(bsr)),y(1..z)),(dp(n),dp(m)),
296      // contains poly check s.t. LT(p) is of the form check(LT(f1),...,LT(fr))
297      def s1 = L[2];
298      map psi = s1,maxideal(1),dom;
299      poly re = p - psi(check);
300      // divide by the maximal power of #[1]
301      if ( (size(#)>0) && (typeof(#[1])=="poly") )
302      {
303        while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
304        {
305          re=re/#[1];
306        }
307      }
308      return(re);
309    }
310    return(p);
311  }
312  //======================2end variant of algorithm=========================
313  //It uses two different commands for elimaination.
314  //if(choose==1):"elimainate"command.
315  //if (choose==2):"nselect" command.
316  else
317  {
318    poly v=product(maxideal(1));
319
320    //------------- change the basering bsr to bsr[@(0),...,@(z)] ----------
321    execute("ring s=("+charstr(basering)+"),("+varstr(basering)+",@(0..z)),dp;");
322    // Ev hier die Reihenfolge der Vars aendern. Dazu muss unten aber entsprechend
323    // geaendert werden:
324    //  execute("ring s="+charstr(basering)+",(@(0..z),"+varstr(basering)+"),dp;");
325
326    //constructs the leading ideal of dom=(p-@(0),dom[1]-@(1),...,dom[z]-@(z))
327    ideal dom=imap(bsr,dom);
328    for (i=1;i<=z;i++)
329    {
331    }
333
334    //---------- eliminate the variables of the basering bsr --------------
335    //i.e. computes dom intersected with K[@(0),...,@(z)].
336
337    if(choose==1)
338    {
339      ideal kern=eliminate(dom,imap(bsr,v));//eliminate does not need a
340                                            //standard basis as input.
341    }
342    if(choose==2)
343    {
344      ideal kern= nselect(groebner(dom),1,n);//"nselect" is combinatorial command
345                                         //which uses the internal command
346                                         // "simplify"
347    }
348
349    //---------  test wether @(0)-h(@(1),...,@(z)) is in ker ---------------
350    // for some poly h and divide by maximal power of q=#[1]
351    poly h;
352    z=size(kern);
353    for (i=1;i<=z;i++)
354    {
355      h=kern[i]/@(0);
356      if (deg(h)==0)
357      {
358        h=(1/h)*kern[i];
359        // define the map psi : s ---> bsr defined by @(i) ---> p,dom[i]
360        setring bsr;
361        map psi=s,maxideal(1),p,dom;
362        poly re=psi(h);
363        // divide by the maximal power of #[1]
364        if ((size(#)>0) && (typeof(#[1])== "poly") )
365        {
366          while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
367          {
368            re=re/#[1];
369          }
370        }
371        return(re);
372      }
373    }
374    setring bsr;
375    return(p);
376  }
377}
378example
379{"EXAMPLE:"; echo = 2;
380  ring r= 0,(x,y),dp;
381  ideal dom =x2,y2,xy-y;
382  poly p=x4+x3y+xy2-y2;
383  sagbiReduction(p,dom);
384  sagbiReduction(p,dom,1);
385  sagbiReduction(p,dom,2);
386}
387
388///////////////////////////////////////////////////////////////////////////////
389static proc completeReduction(poly p,ideal dom,list#)//reduction
390{
391  poly p1=p;
392  poly p2=sagbiReduction(p,dom,#);
393  while (p1!=p2)
394  {
395    p1=p2;
396    p2=sagbiReduction(p1,dom,#);
397  }
398  return(p2);
399}
400///////////////////////////////////////////////////////////////////////////////
401
402static proc completeReduction1(poly p,ideal dom,list #) //tail reduction
403{
404  poly p1,p2,re;
405  p1=p;
406  while(p1!=0)
407  {
408    p2=sagbiReduction(p1,dom,#);
409    if(p2!=p1)
410    {
411      p1=p2;
412    }
413    else
414    {
417    }
418  }
419  return(re);
420}
421
422
423
424///////////////////////////////////////////////////////////////////////////////
425
426proc sagbiNF(id,ideal dom,int k,list#)
427"USAGE: sagbiNF(id,dom,k[,n]); id either poly or ideal,dom ideal, k and n positive intergers.
428RETURN: depends On the type of id; ideal or polynomial.
429@format
430    The integer k determines what kind of s-reduction is performed:
431    - if (k=0) no tail s-reduction is performed.
432    - if (k=1) tail s-reduction is performed.
433    Three Algorthim  variants  are used to perform Subalgebra reduction.
434    The positive integer n determine which variant should be used.
435    n may take the values (0 or default),1 or 2.
436@end format
437NOTE: computation of Subalgebras normal forms may be performed either
438      in polynomial rings or quotient polynomial rings
439EXAMPLE: example sagbiNF; show example "
440{
441  int z;
442  ideal Red;
443  poly re;
444  if(typeof(id)=="ideal")
445  {
446    int i=ncols(id);
447    for(z=1;z<=i;z++)
448    {
449      if(k==0)
450      {
451        id[z]=completeReduction(id[z],dom,#);
452      }
453      else
454      {
455        id[z]=completeReduction1(id[z],dom,#);//tail reduction.
456      }
457    }
458    Red=simplify(id,7);
459    return(Red);
460  }
461  if(typeof(id)=="poly")
462  {
463    if(k==0)
464    {
465      re=completeReduction(id,dom,#);
466    }
467    else
468    {
469      re=completeReduction1(id,dom,#);
470    }
471    return(re);
472  }
473}
474example
475{"EXAMPLE:"; echo = 2;
476 ring r=0,(x,y),dp;
477 ideal I= x2-xy;
478 qring Q=std(I);
479 ideal dom =x2,x2y+y,x3y2;
480 poly p=x4+x2y+y;
481 sagbiNF(p,dom,0);
482 sagbiNF(p,dom,1);// tail subalgebra reduction is performed
483}
484
485
486///////////////////////////////////////////////////////////////////////////////
487
488static proc intRed(id,int k, list #)
489{
490  int i,z;
491  ideal Rest,intRed;
492  z=ncols(id);
493  for(i=1;i<=z;i++)
494  {
495    Rest=id;
496    Rest[i]=0;
497    Rest=simplify(Rest,2);
498    if(k==0)
499    {
500      intRed[i]=completeReduction(id[i],Rest,#);
501    }
502    else
503    {
504      intRed[i]=completeReduction1(id[i],Rest,#);
505    }
506  }
507  intRed=simplify(intRed,7);//1+2+4 in simplify command
508  return(intRed);
509}
510//////////////////////////////////////////////////////////////////////////////
511
512proc sagbi(id,int k,list#)
513"USAGE: sagbi(id,k[,n]); id ideal, k and n positive integers.
514RETURN: A SAGBI basis for the subalgebra defined by the generators of id.
515@format
516    k determine what kind of s-reduction is performed:
517     - if (k=0) no tail s-reduction is performed.
518     - if (k=1) tail s-reduction is performed, and S-interreduced SAGBI basis
519       is returned.
520    Three Algorithm variants are used to perform Subalgebra reduction.
521    The positive interger n determine which variant should be used.
522    n may take the values (0 or default),1 or 2.
523@end format
524NOTE: SAGBI bases computations may be performed either
525      in polynomial rings or quotient polynomial rings.
526EXAMPLE: example sagbi; show example "
527{
528  degBound=0;
529  ideal S,oldS,Red;
530  list L;
531  S=intRed(id,k,#);
532  while(size(S)!=size(oldS))
533  {
534    L=Spoly1(L,S,Red,2);
535    Red=L[1];
536    Red=sagbiNF(Red,S,k,#);
537    oldS=S;
538    S=S+Red;
539  }
540  return(S);
541}
542example
543{ "EXAMPLE:"; echo = 2;
544 ring r= 0,(x,y),dp;
545 ideal I=x2,y2,xy+y;
546 sagbi(I,1,1);
547}
548///////////////////////////////////////////////////////////////////////////////
549proc sagbiPart(id,int k,int c,list #)
550"USAGE: sagbi(id,k,c[,n]); id ideal, k, c and n positive integer.
551RETURN: A partial SAGBI basis for the subalgebra defined by the generators of id.
552@format
553     should stop. k determine what kind of s-reduction is performed:
554     - if (k=0) no tail s-reduction is performed.
555     - if (k=1) tail s-reduction is performed, and S-intereduced SAGBI basis
556      is returned.
557     c  determines, after which turn Sagbi basis computations should stop
558     Three Algorithm variants are used to perform Subalgebra reduction.
559     The positive integer n determines which variant should be used.
560     n may take the values (0 or default),1 or 2.
561@end format
562NOTE:- SAGBI bases computations may be performed either
563       in polynomial rings or quotient polynomial rings.
564     - This version of sagbi procedure is interesting in the case of an Subalgebras
565       with infinte SAGBI basis. In this case, by means of this procedure,
566       we may check for example, if the elements of this basis have a particular form.
567EXAMPLE: example sagbiPart; show example "
568{
569  degBound=0;
570  ideal S,oldS,Red;
571  int counter;
572  list L;
573  S=intRed(id,k,#);
574  while((size(S)!=size(oldS))&&(counter<=c))
575  {
576    L=Spoly1(L,S,Red,2);
577    Red=L[1];
578    Red=sagbiNF(Red,S,k,#);
579    oldS=S;
580    S=S+Red;
581    counter=counter+1;
582  }
583  return(S);
584}
585example
586{ "EXAMPLE:"; echo = 2;
587 ring r= 0,(x,y),dp;
588 ideal I=x,xy-y2,xy2;//the corresponding Subalgebra has an infinte SAGBI basis
589 sagbiPart(I,1,3);// computations should stop after 3 turns.
590}
591//////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.