source: git/Singular/LIB/teachstd.lib @ 2ab830

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