source: git/Singular/LIB/teachstd.lib @ 1a3911

spielwiese
Last change on this file since 1a3911 was 1a3911, checked in by Frank Seelisch <seelisch@…>, 15 years ago
removed some docu errors prior to release 3-1-0 git-svn-id: file:///usr/local/Singular/svn/trunk@11626 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 23.4 KB
Line 
1// $Id:
2//GMG, last modified 28.9.01
3///////////////////////////////////////////////////////////////////////////////
4version="$Id: teachstd.lib,v 1.11 2009-04-06 12:39:02 seelisch Exp $";
5category="Teaching";
6info="
7LIBRARY:  teachstd.lib   Procedures for teaching standard bases
8AUTHOR:                  G.-M. Greuel, greuel@mathematik.uni-kl.de
9
10NOTE: The library is intended to be used for teaching purposes, but not
11      for serious computations. Sufficiently high printlevel allows to
12      control each step, thus illustrating the algorithms at work.
13      The procedures are implemented exactly as described in the book
14      'A SINGULAR Introduction to Commutative Algebra' by G.-M. Greuel and
15      G. Pfister (Springer 2002).
16
17PROCEDURES:
18 ecart(f);              ecart of f
19 tail(f);               tail of f
20 sameComponent(f,g);    test for same module component of lead(f) and lead(g)
21 leadmonomial(f);       leading monomial as poly (also for vectors)
22 monomialLcm(m,n);      lcm of monomials m and n as poly (also for vectors)
23 spoly(f[,1]);          s-polynomial of f [symmetric form]
24 minEcart(T,h);         element g from T of minimal ecart s.t. LM(g)|LM(h)
25 NFMora(i);             normal form of i w.r.t Mora algorithm
26 prodcrit(f,g[,o]);     test for product criterion
27 chaincrit(f,g,h);      test for chain criterion
28 pairset(G);            pairs form G neither satifying prodcrit nor chaincrit
29 updatePairs(P,S,h);    pairset P enlarded by not useless pairs (h,f), f in S
30 standard(id);          standard basis of ideal/module
31 localstd(id);          local standard basis of id using Lazard's method
32
33              [parameters in square brackets are optional]
34";
35
36LIB "poly.lib";
37
38///////////////////////////////////////////////////////////////////////////////
39proc ecart(f)
40"USAGE:   ecart(f); f poly or vector
41RETURN:  the ecart e of f of type int
42EXAMPLE: example ecart; shows an example
43"
44{
45  int e = maxdeg1(f)-maxdeg1(lead(f));
46  return(e);
47}
48example
49{ "EXAMPLE:"; echo = 2;
50   ring r=0,(x,y,z),ls;
51   ecart((y+z+x+xyz)**2);
52   ring s=0,(x,y,z),dp;
53   ecart((y+z+x+xyz)**2);
54}
55///////////////////////////////////////////////////////////////////////////////
56proc leadmonomial(f)
57"USAGE:   leadmonomial(f); f poly or vector
58RETURN:  the leading monomial of f of type poly
59NOTE:    if f is of type poly, leadmonomial(f)=leadmonom(f), if f is of type
60         vector and if leadmonom(f)=m*gen(i) then leadmonomial(f)=m
61EXAMPLE: example leadmonomial; shows an example
62"
63{
64  int e;
65  poly m;
66  if(typeof(f) == "vector")
67  {
68    e=leadexp(f)[nvars(basering)+1];
69    m=leadmonom(f)[e,1];
70  }
71  if(typeof(f) == "poly")
72  {
73    m=leadmonom(f);
74  }
75  return(m);
76}
77example
78{ "EXAMPLE:"; echo = 2;
79   ring s=0,(x,y,z),(c,dp);
80   leadmonomial((y+z+x+xyz)^2);
81   leadmonomial([(y+z+x+xyz)^2,xyz5]);
82}
83///////////////////////////////////////////////////////////////////////////////
84proc tail(f)
85"USAGE:   tail(f); f poly or vector
86RETURN:  f-lead(f), the tail of f of type poly
87EXAMPLE: example tail; shows an example
88"
89{
90  def t = f-lead(f);
91  return(t);
92}
93example
94{ "EXAMPLE:"; echo = 2;
95   ring r=0,(x,y,z),ls;
96   tail((y+z+x+xyz)**2);
97   ring s=0,(x,y,z),dp;
98   tail((y+z+x+xyz)**2);
99}
100///////////////////////////////////////////////////////////////////////////////
101proc sameComponent(f,g)
102"USAGE:   sameComponent(f,g); f,g poly or vector
103RETURN:  1 if f and g are of type poly or if f and g are of type vector and
104         their leading monomials involve the same module component,
105         0 if not
106EXAMPLE: example sameComponent; shows an example
107"
108{
109  if(typeof(f) != typeof(g))
110  {
111    ERROR("** arguments must be of same type");
112  }
113  if(typeof(f) == "vector")
114  {
115     if( leadexp(f)[nvars(basering)+1] !=  leadexp(g)[nvars(basering)+1] )
116     {
117       return(0);
118     }
119  }
120  return(1);
121}
122example
123{ "EXAMPLE:"; echo = 2;
124   ring r=0,(x,y,z),dp;
125   sameComponent([y+z+x,xyz],[z2,xyz]);
126   sameComponent([y+z+x,xyz],[z4,xyz]);
127   sameComponent(y+z+x+xyz, xy+z5);
128}
129///////////////////////////////////////////////////////////////////////////////
130proc monomialLcm(m,n)
131"USAGE:   monomialLcm(m,n); m,n of type poly or vector
132RETURN:  least common multiple of leading monomials of m and n, of type poly
133NOTE:    if m = (x1...xr)^(a1,...,ar)*gen(i) (gen(i)=1 if m is of type poly)
134         and n = (x1...xr)^(b1,...,br)*gen(j), then the proc returns
135         (x1,...,xr)^(max(a1,b1),...,max(ar,br)) if i=j and 0 if i!=j.
136EXAMPLE: example monomialLcm; shows an example
137"
138{
139  if(typeof(n) != typeof(m))
140  {
141    ERROR("** arguments must be of same type");
142  }
143
144  poly u ;
145  if(sameComponent(m,n) == 0)   //leading term of vectors involve
146  {                             //different module components
147     return(u);
148  }
149
150  intvec v = leadexp(m);        //now start to compute lcm
151  intvec w = leadexp(n);
152  u=1;
153  int i;
154  for (i=1; i<=nvars(basering); i++)
155  {
156      if(v[i]>=w[i])
157      {
158         u = u*var(i)**v[i];
159      }
160      else
161      {
162         u = u*var(i)**w[i];
163      }
164   }
165   return(u);
166}
167example
168{ "EXAMPLE:"; echo = 2;
169   ring r=0,(x,y,z),ds;
170   monomialLcm(xy2,yz3);
171   monomialLcm([xy2,xz],[yz3]);
172   monomialLcm([xy2,xz3],[yz3]);
173}
174///////////////////////////////////////////////////////////////////////////////
175proc spoly(f,g,list #)
176"USAGE:   spoly(f,g[,s]); f,g poly or vector, s int
177RETURN:  the s-polynomial of f and g, of type poly or vector
178         if s!=0 the symmetric s-polynomial (without division) is returned
179EXAMPLE: example spoly; shows an example
180"
181{
182  if(typeof(f) != typeof(g))
183  {
184    ERROR("** arguments must be of same type");
185  }
186  if(size(#) == 0)
187  {
188     #[1]=0;
189  }
190
191  int e;
192  poly o = monomialLcm(f,g);
193
194  if( o == 0)                    //can only happen, if vectors f and g involve
195  {                              //different module components
196    vector sp;
197    return(sp);
198  }
199
200  poly m=(o/leadmonomial(f));         //compute the leading monomial as poly
201  poly n=(o/leadmonomial(g));
202
203  f = m * f;
204  g = n * g; // now they have the same LM!
205
206  if (#[1]==0)                    //the asymmetric s-poly
207  {
208    def sp = f - (leadcoef(f)/leadcoef(g))*g;
209  }
210  else                            //the symmetric s-poly, avoiding division
211  {
212    def sp = leadcoef(g)*f - leadcoef(f)*g;
213  }
214  return(sp);
215}
216example
217{ "EXAMPLE:"; echo = 2;
218   ring r=0,(x,y,z),ls;
219   spoly(2x2+x2y,3y3+xyz);
220   ring s=0,(x,y,z),(c,dp);
221   spoly(2x2+x2y,3y3+xyz);
222   spoly(2x2+x2y,3y3+xyz,1);             //symmetric s-poly without division
223   spoly([5x2+x2y,z5],[x2,y3,y4]);       //s-poly for vectors
224}
225///////////////////////////////////////////////////////////////////////////////
226proc minEcart(T,h)
227"USAGE:   minEcart(T,h); T ideal or module, h poly or vector
228RETURN:  element g from T such that leadmonom(g) divides leadmonom(h)@*
229         ecart(g) is minimal with this property (if T != 0);
230         return 0 if T is 0 or h = 0
231EXAMPLE: example minEcart; shows an example
232"
233{
234  int i,k,e;
235  if (size(T)==0 or h==0 )       //trivial cases
236  {
237     h = 0;
238     return(h);
239  }
240//---- check whether some element in T has the same module component as h ----
241  int v;
242  intvec w;
243  T = simplify(T,2);
244
245  if (typeof(h) == "vector")
246  {
247     e=1;
248     v = leadexp(h)[nvars(basering)+1];
249     for( i=1; i<=size(T); i++)
250     {
251       w[i]=leadexp(T[i])[nvars(basering)+1];
252       if(v == w[i])
253       {
254         e=0;         //some element in T involves the same component as h
255       }
256     }
257  }
258  if ( e == 1 )       //no element in T involves the same component as h
259  {
260     h = 0;
261     return(h);
262  }
263
264  if (typeof(h) == "poly")   //for polys v=w[i] for all i
265  {
266    v = 1;
267    w[size(T)]=0;
268    w=w+1;
269  }
270
271//------ check whether for some g in T leadmonom(g) divides leadmonom(h) -----
272  def g = T[1];
273  poly f = leadmonomial(h);
274
275  for( i=1; i<=size(T); i++)
276  {
277    if( f/leadmonomial(T[i]) != 0 and v==w[i] )
278    {
279       g=T[i];
280       k=i;
281       break;
282    }
283  }
284  if (k == 0)        //no leadmonom(g) divides leadmonom(h)
285  {
286     g = 0;
287     return(g);
288  }
289//--look for T[i] with minimal ecart s.t.leadmonom(T[i]) divides leadmonom(h)--
290
291  for( i=k+1; i<=size(T); i++)
292  {
293    if( f/leadmonomial(T[i]) != 0 and v==w[i] )
294    {
295       if (ecart(T[i]) < ecart(g))
296       {
297         g=T[i];
298        }
299     }
300  }
301  return(g);
302}
303example
304{ "EXAMPLE:"; echo = 2;
305   ring R=0,(x,y,z),dp;
306   ideal T = x2y+x2,y3+xyz,xyz2+z4;
307   poly h = x2y2z2+x5+yx3+z6;
308   minEcart(T,h);"";
309   ring S=0,(x,y,z),(c,ds);
310   module T = [x2+x2y,y2],[y3+xyz,x3-z3],[x3y+z4,0,x2];
311   vector h = [x3y+x5+x2y2z2+z6,x3];
312   minEcart(T,h);
313}
314///////////////////////////////////////////////////////////////////////////////
315proc NFMora(f,G,list #)
316"USAGE:   NFMora(f,G[,s]); f poly or vector,G ideal or module, s int
317RETURN:  the Mora normal form of f w.r.t. G, same type as f
318         if s!=0 the symmetric s-polynomial (without division) is used
319NOTE:    Show comments if printlevel > 0, pauses computation if printlevel > 1
320EXAMPLE: example NFMora; shows an example
321"
322{
323  if(size(#) == 0)
324  {
325     #[1]=0;
326  }
327
328  int y = printlevel - voice + 2;
329  if( f==0 or size(G) ==0 )
330  {
331     if (y>0)
332     {
333        "// 1st or 2nd argument 0";
334     }
335     return(f);
336  }
337
338  int i,e;
339  def h = f;
340  def T = G;
341// -------------------- start with f to be reduced by G --------------------
342  if (y>0)
343  {"";
344   "// Input for NFMora is (f,T):";
345   "// f:"; f;
346   "// T:"; T;
347   "// Reduce f with T, eventually enlarging T for local ordering";
348  }
349
350// ----------------------- enter the reduction loop ------------------------
351  def g = minEcart(T,h);
352  while (h!=0 and g!=0)
353  {
354     if (y>0)
355     {  "";
356        "//                Reduction-step in NFMora:",i;
357        "// h = (f after",i,"reductions) reduction with g from T:";
358        "// g = element of minimal ecart in T s.t. LM(g)|LM(h):";
359        "// h:";h;
360        "// g:";g;
361     }
362     if (y>1)
363     {
364        pause("press <return> to continue");
365        "// T, set used for reduction:"; T;
366        pause("press <return> to continue");
367     }
368     e=0;
369     if( ecart(g) > ecart(h) )
370     {
371        T=T,h;  e=1;
372     }
373     if (y>0 )
374     {
375        "// T-set enlarged for next reduction? (yes/no = 1/0):  ", e;
376        if( e==1 )
377        {
378          "// T-set for next reduction got enlarged by h:";
379          "// h:";h;
380          if (y>1)
381          {
382             pause("press <return> to continue");
383          }
384        }
385     }
386     h = spoly(h,g,#[1]);
387     g = minEcart(T,h);
388     i=i+1;
389  }
390  if(y>0)
391  { "";
392    "// normal form is:";
393  }
394  return(h);
395}
396example
397{ "EXAMPLE:"; echo = 2;
398   ring r=0,(x,y,z),dp;
399   poly f = x2y2z2+x5+yx3+z6-3y3;
400   ideal G = x2y+x2,y3+xyz,xyz2+z6;
401   NFMora(f,G);"";
402   ring s=0,(x,y,z),ds;
403   poly f = x3y+x5+x2y2z2+z6;
404   ideal G = x2+x2y,y3+xyz,x3y2+z4;
405   NFMora(f,G);"";
406   vector v = [f,x2+x2y];
407   module M = [x2+x2y,f],[y3+xyz,y3],[x3y2+z4,z2];
408   NFMora(v,M);
409}
410///////////////////////////////////////////////////////////////////////////////
411proc prodcrit(f,g,list #)
412"USAGE:   prodcrit(f,g[,o]); f,g poly or vector, and optional int argument o
413RETURN:  1 if product criterion applies in the same module component,
414         2 if lead(f) and lead(g) involve different components, 0 else
415NOTE:    if product criterion applies we can delete (f,g) from pairset.
416         This procedure returns 0 if o is given and is a positive integer, or
417         you may set the attribute \"default_arg\" for prodcrit to 1.
418EXAMPLE: example prodcrit; shows an example
419"
420{
421// ------------------ check for optional disabling argument -------------
422  if( size(#) > 0 )
423  {// "size(#): ", size(#);   "typeof(#[1]): ", typeof(#[1]);
424   
425    if( typeof(#[1]) == "int" )
426    {// "#[1] = int ", #[1];
427      if( #[1] > 0 )
428      {       
429        return(0);
430      }
431    }
432  }
433
434// ------------------- product criterion for polynomials ---------------------
435  if(typeof(f)=="poly")
436  {
437    if( monomialLcm(f,g)==leadmonom(f)*leadmonom(g))
438    {
439      return(1);
440    }
441    return(0);
442  }
443// ------------------- product criterion for modules ---------------------
444  if(sameComponent(f,g)==1)
445  {
446     if( monomialLcm(f,g)==leadmonomial(f)*leadmonomial(g) )
447     {
448         int c = leadexp(f)[nvars(basering)+1];  //component involving lead(f)
449         if((f-f[c]*gen(c))-(g-g[c]*gen(c))==0)  //other components are 0
450         {
451           return(1);
452         }
453     }
454     return(0);
455  }
456  return(2);
457}
458example
459{ "EXAMPLE:"; echo = 2;
460   ring r=0,(x,y,z),dp;
461   poly f = y3z3+x5+yx3+z6;
462   poly g = x5+yx3;
463   prodcrit(f,g);
464   vector v = x3z2*gen(1)+x3y*gen(1)+x2y*gen(2);
465   vector w = y4*gen(1)+y3*gen(2)+xyz*gen(1);
466   prodcrit(v,w);
467}
468///////////////////////////////////////////////////////////////////////////////
469proc chaincrit(f,g,h)
470"USAGE:   chaincrit(f,g,h); f,g,h poly or module
471RETURN:  1 if chain criterion applies, 0 else
472NOTE:    if chain criterion applies to f,g,h we can delete (g,h) from pairset
473EXAMPLE: example chaincrit; shows an example
474"
475{
476  if(sameComponent(f,g) and sameComponent(f,h))
477  {
478    if( monomialLcm(g,h)/leadmonomial(f) !=0 )
479    {
480      return(1);
481    }
482  }
483  return(0);
484}
485example
486{ "EXAMPLE:"; echo = 2;
487   ring r=0,(x,y,z),dp;
488   poly f = x2y2z2+x5+yx3+z6;
489   poly g = x5+yx3;
490   poly h = y2z5+x5+yx3;
491   chaincrit(f,g,h);
492   vector u = [x2y3-z2,x2y];
493   vector v = [x2y2+z2,x2-y2];
494   vector w = [x2y4+z3,x2+y2];
495   chaincrit(u,v,w);
496}
497///////////////////////////////////////////////////////////////////////////////
498proc pairset(G)
499"USAGE:   pairset(G); G ideal or module
500RETURN:  list L,
501         L[1] = the pairset of G as list (not containing pairs for
502         which the product or the chain criterion applies),
503         L[2] = intvec v, v[1]= # product criterion, v[2]= # chain criterion
504EXAMPLE: example pairset; shows an example
505"
506{
507   int i,j,k,s,c,ccrit,pcrit,pr;
508   int y = printlevel - voice + 2;
509   G = simplify(G,10);
510   def g = G;
511   ideal pair;
512   list P,I;      //P=pairlist of G, I=list of corresponding indices of pairs
513   for (i=1; i<=size(G); i++)
514   {
515      for(j = i+1; j<=size(G); j++)
516      {
517         pr = prodcrit(G[i],G[j]);       //check first product criterion
518         if( pr != 0 )
519         {
520            pcrit=pcrit+(pr==1);
521         }
522         else
523         {
524            s = size(P);                //now check chain criterion
525            for(k=1; k<=s; k++)
526            {
527              if( i==I[k][2] )
528              {
529                if ( chaincrit(P[k][1],P[k][2],G[j]) )
530                {                     //need not include (G[i],G[j]) in P
531                  c=1; ccrit=ccrit+1;
532                  break;
533                }
534              }
535              if( j==I[k][1] and c==0 )
536              {
537                "########### enter pairset2 #############";
538                if ( chaincrit(G[i],P[k][1],P[k][2]) )
539                {                     //can delete P[k]=(P[k][1],P[k][2])
540                  ccrit=ccrit+1;
541                  P = delete(P,k);
542                  s = s-1;
543                }
544              }
545            }
546            if ( c==0 )
547            {
548               g = G[i],G[j];
549               P[s+1]=g;
550               I[s+1]=intvec(i,j);
551            }
552            c=0;
553         }
554      }
555   }
556   if (y>0)
557   {  "";
558      "// product criterion:",pcrit," chain criterion:",ccrit;
559   }
560   intvec v = pcrit,ccrit;
561   P=P,v;
562   return(P);
563}
564example
565{ "EXAMPLE:"; echo = 2;
566   ring r=0,(x,y,z),dp;
567   ideal G = x2y+x2,y3+xyz,xyz2+z4;
568   pairset(G);"";
569   module T = [x2y3-z2,x2y],[x2y2+z2,x2-y2],[x2y4+z3,x2+y2];
570   pairset(T);
571}
572///////////////////////////////////////////////////////////////////////////////
573proc updatePairs(P,S,h)
574"USAGE:   updatePairs(P,S,h); P list, S ideal or module, h poly or vector@*
575         P a list of pairs of polys or vectors (obtained from pairset)
576RETURN:  list Q,
577         Q[1] = the pairset P enlarged  by all pairs (f,h), f from S,
578         without pairs for which the product or the chain criterion applies@*
579         Q[2] = intvec v, v[1]= # product criterion, v[2]= # chain criterion
580EXAMPLE: example updatePairs; shows an example
581"
582{
583   int i,j,k,s,r,c,ccrit,pcrit,pr;
584   int y = printlevel - voice + 2;
585   ideal pair;
586   list Q = P;           //Q will become enlarged pairset
587   s = size(P);
588   r = size(Q);          //r will grow with Q
589   list R;
590   def g = S;            //give g the correct type ideal/module
591   for (i=1; i<=size(S); i++)
592   {
593      pr =  prodcrit(h,S[i]);
594      if( pr != 0 )     //check product criterion
595      {
596        pcrit=pcrit+(pr==1);       //count product criterion in same component
597      }
598      else
599      {                 //prodcrit did not apply, check for chaincrit
600         r=size(Q);
601         for(k=1; k<=r; k++)
602         {
603           if( Q[k][2]==S[i] )    //S[i]=Q[k][2]
604           {
605             if( chaincrit(Q[k][1],S[i],h) )
606             {                   //can forget (S[i],h)
607                c=1; ccrit=ccrit+1;
608                break;
609             }
610           }
611         }
612         if ( c==0 )
613         {
614            g = S[i],h;        //add pair (S[i],h)
615            Q[r+1] = g;
616         }
617         c=0;
618      }
619   }
620   if (y>0)
621   {  "";;
622      "// product criterion:",pcrit," chain criterion:",ccrit;
623   }
624   intvec v = pcrit,ccrit;
625   Q = Q,v;
626   return(Q);
627}
628example
629{ "EXAMPLE:"; echo = 2;
630   ring R1=0,(x,y,z),(c,dp);
631   ideal S = x2y+x2,y3+xyz;
632   poly h = x2y+xyz;
633   list P = pairset(S)[1];
634   P;"";
635   updatePairs(P,S,h);"";
636   module T = [x2y3-z2,x2y],[x2y4+z3,x2+y2];
637   P = pairset(T)[1];
638   P;"";
639   updatePairs(P,T,[x2+x2y,y3+xyz]);
640}
641///////////////////////////////////////////////////////////////////////////////
642proc standard(id, list #)
643"USAGE:   standard(i[,s]); id ideal or module, s int
644RETURN:  a standard basis of id, using generalized Mora's algorithm
645         which is Buchberger's algorithm for global monomial orderings.
646         If s!=0 the symmetric s-polynomial (without division) is used
647NOTE:    Show comments if printlevel > 0, pauses computation if printlevel > 1
648EXAMPLE: example standard; shows an example
649"
650{
651  if(size(#) == 0)
652  {
653     #[1]=0;
654  }
655
656   def S = id;                 //S will become the standard basis of id
657   def h = S[1];
658   int i,z;
659   int y = printlevel - voice + 2;
660   if(y>0)
661   {  "";
662      "// the set S, to become a standard basis:"; S;
663      if(y>1)
664      {
665         "// create pairset, i.e. pairs from S,";
666         "// after application of product and chain criterion";
667      }
668   }
669   list P = pairset(S);        //create pairset of S=id
670   intvec v = P[2];
671   P = P[1];
672//-------------------------- Main loop in SB lgorithm ----------------------
673   while (size(P) !=0)
674   {  z=z+1;
675      if(y>0)
676      {  "";
677         "//              Enter NFMora for next pair, count",z;
678         "// size of partial standard basis S:    (",size(S),")";
679         "// number of pairs of S after updating: (",size(P),")";
680         if(y>1)
681         {
682            "// 1st pair of new pairset:"; P[1];
683            "// set T=S used for reduction:";S;
684            "// apply NFMora to (spoly,S), spoly = spoly(1st pair)";
685         }
686      }
687      //-------------------- apply NFMora = Mora's normal form -------------
688      h = spoly(P[1][1],P[1][2],#[1]);
689      if(y>1)
690      {
691         "// spoly:";h;
692      }
693      h = NFMora(h,S,#[1]);
694      if(h==0)           //normal form is 0
695      {
696         if(y==1)
697         {
698            "// pair has reduced to 0";
699         }
700          if(y>1)
701         { h;"";
702           pause("press <return> to continue");
703         }
704     }
705      P = delete(P,1);   //spoly of pair reduced to 0, pair can be deleted
706      //--- spoly of pair did not reduce to 0, update S and paiset of S ----
707      if( h != 0)
708      {
709         if(y==1)
710         {
711            "// ** new spoly in degree **:", deg(h);
712         }
713         if(y>1)
714         { h;"";
715           pause("press <return> to continue");
716           "// update pairset";
717         }
718        P=updatePairs(P,S,h);    //update P (=paisert of S)
719        v=v+P[2];                //with useful pairs (g,h), g from S
720         P=P[1];
721         S=S,h;                  //update S, will become the standard basis
722      }
723   }
724//------------------------------ finished ---------------------------------
725   if( find(option(),"prot") or y>0 )
726   {  "";                //show how often prodcrit and chaincrit applied
727      "// product criterion:",v[1]," chain criterion:",v[2];
728      "";
729      "// Final standard basis:";
730   }
731   return(S);
732}
733example
734{ "EXAMPLE:"; echo = 2;
735   ring r=0,(x,y,z),dp;
736   ideal G = x2y+x2,y3+xyz,xyz2+z4;
737   standard(G);"";
738   ring s=0,(x,y,z),(c,ds);
739   ideal G = 2x2+x2y,y3+xyz,3x3y+z4;
740   standard(G);"";
741   standard(G,1);"";           //use symmetric s-poly without division
742   module M = [2x2,x3y+z4],[3y3+xyz,y3],[5z4,z2];
743   standard(M);
744}
745///////////////////////////////////////////////////////////////////////////////
746proc localstd (id)
747"USAGE:   localstd(id); id = ideal
748RETURN:  A standard basis for a local degree ordering, using Lazard's method.
749NOTE:    The procedure homogenizes id w.r.t. a new 1st variable local@t,
750         computes a SB w.r.t. (dp(1),dp) and substitutes local@t by 1.
751         Hence the result is a SB with respect to an ordering which sorts
752         first w.r.t. the subdegree of the original variables and then refines
753         it with dp. This is the local degree ordering ds.
754         localstd may be used in order to avoid cancellation of units and thus
755         to be able to use option(contentSB) also for local orderings.
756EXAMPLE: example localstd; shows an example
757"
758{
759  int ii;
760  def bas = basering;
761  execute("ring  @r_locstd
762     =("+charstr(bas)+"),(local@t,"+varstr(bas)+"),(dp(1),dp);");
763  ideal id = imap(bas,id);
764  ideal hid = homog(id,local@t);
765  hid = std(hid);
766  hid = subst(hid,local@t,1);
767  setring bas;
768  def hid = imap(@r_locstd,hid);
769  attrib(hid,"isSB",1);
770  kill @r_locstd;
771  return(hid);
772}
773example
774{ "EXAMPLE:"; echo = 2;
775   ring R = 0,(x,y,z),ds;
776   ideal i  = xyz+z5,2x2+y3+z7,3z5+y5;
777   localstd(i);
778}
779///////////////////////////////////////////////////////////////////////////////
780
781/*
782// some examples:
783LIB"teachstd.lib";
784option(prot); printlevel=3;
785ring r0 = 0,(t,x,y),lp;
786ideal i = x-t2,y-t3;
787standard(i);
788
789printlevel=1;
790standard(i);
791
792option(prot); printlevel =1;
793ring r1 = (0,a,b),(x,y,z),(c,ds);
794module M = [ax2,bx3y+z4],[a3y3+xyz,by3],[5az4,(a+b)*z2];
795module N1= std(M);
796module N2 = standard(M,1);
797NF(lead(N2),lead(N1));
798NF(lead(N1),lead(N2));rom T
799ring r2 = 0,(x,y,z),dp;
800ideal I = x2y+x2,y3+xyz,xyz2+z4;
801option(prot);
802int tt = timer;
803ideal J = standard(I);
804timer -tt;                 //4sec, product criterion: 9  chain criterion: 6
805ideal J1 = std(I);
806NF(lead(J),lead(J1));
807NF(lead(J1),lead(J));
808
809ring r3 = 0,(x,y,z),ds;
810poly f = x3*y4+z5+xyz;
811ideal I = f,jacob(f);
812option(prot);
813int tt = timer;
814ideal J = standard(I);
815timer -tt;                 //3sec, product criterion: 1  chain criterion: 3
816ideal J1 = std(I);
817NF(lead(J),lead(J1));
818NF(lead(J1),lead(J));
819
820//Becker:
821ring r4 = 32003,(x,y,z),lp;
822ideal I = x3-1, y3-1,
823-27x3-243x2y+27x2z-729xy2+162xyz-9xz2-729y3+243y2z-27yz2+z3-27;
824option(prot);
825int tt = timer;
826ideal J = standard(I);
827timer -tt;                 //201sec, product criterion: 42  chain criterion: 33
828ideal J1 = std(I);
829NF(lead(J),lead(J1));
830NF(lead(J1),lead(J));
831
832//Alex
833ring r5 = 32003,(x,y,z,t),dp;
834ideal I =
8352t3xy2z+x2ty+2x2y,
8362tz+y3x2t+z2t3y2x;
837option(prot); printlevel =1;
838ideal J1 = std(I);
839int tt = timer;
840ideal J = standard(I);
841timer -tt;                 //15sec product criterion: 0  chain criterion: 12
842NF(lead(J),lead(J1));
843NF(lead(J1),lead(J));
844
845ring r6 = 32003,(x,y,z,t),dp; //is already SB for ds, for dp too long
846ideal I=
8479x8+y7t3z4+5x4y2t2+2xy2z3t2,
8489y8+7xy6t+2x5y4t2+2x2yz3t2,
8499z8+3x2y3z2t4;
850option(prot);
851int tt = timer;
852ideal J = standard(I);
853timer -tt;                 //0sec, product criterion: 3  chain criterion: 0
854ideal J1 = std(I);
855NF(lead(J),lead(J1));
856NF(lead(J1),lead(J));
857
858
859*/
860
Note: See TracBrowser for help on using the repository browser.