# source:git/Singular/LIB/presolve.lib@35ceb0

spielwiese
Last change on this file since 35ceb0 was 35ceb0, checked in by Hans Schönemann <hannes@…>, 18 years ago
• Property mode set to `100644`
File size: 46.3 KB
Line
1//last change: 13.02.2001 (Eric Westenberger)
2///////////////////////////////////////////////////////////////////////////////
3version="\$Id: presolve.lib,v 1.23 2005-05-04 09:31:39 Singular Exp \$";
4category="Symbolic-numerical solving";
5info="
6LIBRARY:  presolve.lib     Pre-Solving of Polynomial Equations
7AUTHOR:   Gert-Martin Greuel, email: greuel@mathematik.uni-kl.de,
8
9PROCEDURES:
10 degreepart(id,d1,d2);  elements of id of total degree >= d1 and <= d2
11 elimlinearpart(id);    linear part eliminated from id
12 elimpart(id[,n]);      partial elimination of vars [among first n vars]
13 elimpartanyr(i,p);     factors of p partially eliminated from i in any ring
14 fastelim(i,p[..]);     fast elimination of factors of p from i [options]
15 faststd(id,[..]);      fast std for ideal id [options]
17 hilbvec(id[,c,o]);     intvec of Hilberseries of id [in char c and ord o]
18 linearpart(id);        elements of id of total degree <=1
19 tolessvars(id[,]);     maps id to new basering having only vars occuring in id
20 solvelinearpart(id);   reduced std-basis of linear part of id
21 sortandmap(id[..]);   map to new basering with vars sorted w.r.t. complexity
22 sortvars(id[n1,p1..]); sort vars w.r.t. complexity in id [different blocks]
23 valvars(id[..]);       valuation of vars w.r.t. to their complexity in id
24 idealsSimplify(id);    ideal I = eliminate(Id,m) m is a product of variables
25                       which are only linearly involved in the generators of id
26 idealSplit(id,tF,fS);  a list of ideals such that their intersection
27                        has the same radical as id
28                       ( parameters in square brackets [] are optional)
29";
30
31LIB "inout.lib";
32LIB "general.lib";
33LIB "matrix.lib";
34LIB "ring.lib";
35LIB "elim.lib";
36///////////////////////////////////////////////////////////////////////////////
37proc shortid (id,int n,list #)
38"USAGE:   shortid(id,n[,e]); id= ideal/module, n,e=integers
39RETURN:  - if called with two arguments or e=0:
40@*       same type as id, containing generators of id having <= n terms.
41@*       - if called with three arguments and e!=0:
42@*       a list L:
43@*       L[1]: same type as id, containing generators of id having <= n terms.
44@*       L[2]: number of corresponding generator of id
45NOTE:    May be used to compute partial standard basis in case id is to hard
46EXAMPLE: example shortid; shows an example
47"
48{
49  intvec v;
50  int ii;
51  for(ii=1; ii<=ncols(id); ii++)
52  {
53   if (size(id[ii]) <=n and id[ii]!=0 )
54   {
55     v=v,ii;
56   }
57   if (size(id[ii]) > n )
58   {
59       id[ii]=0;
60   }
61  }
62  if( size(v)>1 )
63  {
64    v = v[2..size(v)];
65  }
66  id = simplify(id,2);
67  list L = id,v;
68  if ( size(#)==0 )
69  {
70    return(id);
71  }
72  if ( size(#)!=0 )
73  {
74    if(#[1]==0)
75    {
76      return(id);
77    }
78    if(#[1]!=0)
79    {
80      return(L);
81    }
82  }
83}
84example
85{ "EXAMPLE:"; echo = 2;
86   ring s=0,(x,y,z,w),dp;
87   ideal i = (x3+y2+yw2)^2,(xz+z2)^2,xyz-w2-xzw;
88   shortid(i,3);
89}
90///////////////////////////////////////////////////////////////////////////////
91
92proc degreepart (id,int d1,int d2,list #)
93"USAGE:   degreepart(id,d1,d2[,v]);  id=ideal/module, d1,d1=integers, v=intvec
94RETURN:  generators of id of [v-weighted] total degree >= d1 and <= d2
95         (default: v = 1,...,1)
96EXAMPLE: example degreepart; shows an example
97"
98{
99   def dpart = id;
100   int s,ii = size(id),0;
101   if ( size(#)==0 )
102   {
103      for ( ii=1; ii<=s; ii=ii+1 )
104      {
105         dpart[ii] = (jet(id[ii],d1-1)==0)*(id[ii]==jet(id[ii],d2))*id[ii];
106      }
107   }
108   else
109   {
110      for ( ii=1; ii<=s; ii=ii+1 )
111      {
112         dpart[ii]=(jet(id[ii],d1-1,#[1])==0)*(id[ii]==jet(id[ii],d2,#[1]))*id[ii];
113      }
114   }
115   return(simplify(dpart,2));
116}
117example
118{ "EXAMPLE:"; echo = 2;
119   ring r=0,(x,y,z),dp;
120   ideal i=1+x+x2+x3+x4,3,xz+y3+z8;
121   degreepart(i,0,4);
122   module m=[x,y,z],x*[x3,y2,z],[1,x2,z3,0,1];
123   intvec v=2,3,6;
124   show(degreepart(m,8,8,v));
125}
126///////////////////////////////////////////////////////////////////////////////
127
128proc elimlinearpart (ideal i,list #)
129"USAGE:   elimlinearpart(i[,n]);  i=ideal, n=integer,@*
130          default: n=nvars(basering)
131RETURN:   list L with 5 entries:
132  @format
133  L[1]: (interreduced) ideal obtained from i by substituing
134        from the first n variables those, which appear in a linear part
135        of i, by putting this part into triangular form
136  L[2]: ideal of variables which have been substituted
137  L[3]: ideal, j-th element defines substitution of j-th var in [2]
138  L[4]: ideal of variables of basering, eliminated ones are set to 0
139  L[5]: ideal, describing the map from the basering to itself such that
140        L[1] is the image of i
141  @end format
142NOTE:    the procedure does always interreduce the ideal i internally w.r.t.
143         ordering dp.
144EXAMPLE: example elimlinearpart; shows an example
145"
146{
147   int ii,n,fi,k;
148   string o, newo;
149   intvec getoption = option(get);
150   option(redSB);
151   def P = basering;
152   n = nvars(P);
153//--------------- replace ordering by dp-ordering if necessary ----------------
154   o = "dp("+string(n)+")";
155   fi = find(ordstr(P),o);
156   if( fi == 0 or find(ordstr(P),"a") != 0 )
157   {
158      execute("ring newP = ("+charstr(P)+"),("+varstr(P)+"),dp;");
159      ideal i = imap(P,i);
160   }
161//---------------------------------- start ------------------------------------
162   if ( size(#)!=0 ) {  n=#[1]; }
163   ideal maxi,rest = maxideal(1),0;
164   if ( n < nvars(P) ) { rest = maxi[n+1..nvars(P)]; }
165   attrib(rest,"isSB",1);
166//-------------------- interreduce and find linear part  ----------------------
167// interred does the only real work. Because of ordering dp the linear part is
168// within the first elements after interreduction
169// **: perhaps Bareiss to constant matrix of linear part instead of interred,
170// and/or for big systems, interred only those generators of id
171// which do not contain elements not to be eliminated
172
173   ideal id = interred(i);
174
175   if(id[1] == 1 )                          //special case of unit ideal
176   {
177     setring P;
178     list L = ideal(1), ideal(0), ideal(0), maxideal(1), maxideal(1);
179     option(set,getoption);
180     return(L);
181   }
182
183   for ( ii=1; ii<=size(id); ii++ )
184   {
185      if( deg(id[ii]) > 1) { break; }
186      k=ii;
187   }
188   if( k == 0 )       { ideal lin; }
189   else               { ideal lin = id[1..k];}
190   if( k < size(id) ) { id = id[k+1..size(id)]; }
191   else               { id = 0; }
192//----- remove generators from lin containing vars not to be eliminated  ------
193   if ( n < nvars(P) )
194   {
195      for ( ii=1; ii<=size(lin); ii++ )
196      {
197         if ( reduce(lead(lin[ii]),rest) == 0 )
198         {
199            id=lin[ii],id;
200            lin[ii] = 0;
201         }
202      }
203   }
204   lin = simplify(lin,3);
205   attrib(lin,"isSB",1);
207   attrib(eva,"isSB",1);
208   ideal neva = reduce(maxideal(1),eva);
209//------------------ go back to original ring end return  ---------------------
210   if( fi == 0 or find(ordstr(P),"a") != 0 )
211   {
212      setring P;
213      ideal id = imap(newP,id);
214      ideal eva = imap(newP,eva);
215      ideal lin = imap(newP,lin);
216      ideal neva = imap(newP,neva);
217   }
218
219   eva = eva[ncols(eva)..1];  // sorting according to variables in basering
220   lin = lin[ncols(lin)..1];
221   ideal phi= neva;
222   k = 1;
223   for( ii=1; ii<=n; ii++ )
224   {
225      if( neva[ii] == 0 )
226      {
227         phi[ii] = eva[k]-lin[k];
228         k=k+1;
229      }
230   }
231   list L = id, eva, lin, neva, phi;
232   option(set,getoption);
233   return(L);
234}
235example
236{ "EXAMPLE:"; echo = 2;
237   ring s=0,(x,y,z),dp;
238   ideal i = x3+y2+z,x2y2+z3,y+z+1;
239  elimlinearpart(i);
240}
241///////////////////////////////////////////////////////////////////////////////
242
243proc elimpart (ideal i,list #)
244"USAGE:   elimpart(i [,n,e] );  i=ideal, n,e=integers
245         n   : only the first n vars are considered for substitution,@*
246         e =0: substitute from linear part of i (same as elimlinearpart)@*
247         e!=0: eliminate also by direct substitution@*
248         (default: n = nvars(basering), e = 1)
249RETURN:  list of 5 objects:
250  @format
251  [1]: ideal obtained by substituting from the first n variables those
252       from i, which appear in the linear part of i (or, if e!=0, which
253       can be expressed directly in the remaining vars)
254  [2]: ideal, variables which have been substituted
255  [3]: ideal, i-th element defines substitution of i-th var in [2]
256  [4]: ideal of variables of basering, substituted ones are set to 0
257  [5]: ideal, describing the map from the basering, say k[x(1..m)], to
258       itself onto k[..variables fom [4]..] and [1] is the image of i
259  @end format
260  The ideal i is generated by [1] and [3] in k[x(1..m)], the map [5]
261  maps [3] to 0, hence induces an isomorphism
262  @format
263            k[x(1..m)]/i -> k[..variables fom [4]..]/[1]
264  @end format
265NOTE:    If the basering has ordering (c,dp), this is faster for big ideals,
266         since it avoids internal ring change and mapping.
267EXAMPLE: example elimpart; shows an example
268"
269{
270   def P = basering;
271   int n,e = nvars(P),1;
272   if ( size(#)==1 ) {  n=#[1]; }
273   if ( size(#)==2 ) {  n=#[1]; e=#[2];}
274//----------- interreduce linear part with proc elimlinearpart ----------------
275// lin = ideal after interreduction of linear part
276// eva = eliminated (substituted) variables
277// sub = polynomials defining substitution
278// neva= not eliminated variables
279// phi = map describing substitution
280   list L = elimlinearpart(i,n);
281   ideal lin, eva, sub, neva, phi = L[1], L[2], L[3], L[4], L[5];
282//-------- direct substitution of variables if possible and if e!=0 -----------
283// first find terms lin1 of pure degree 1 in each poly of lin
284// k1 = pure degree 1 part, k1+k2 = those polys of lin which contained a pure
285// degree 1 part.
286// Then go to ring newP with ordering c,dp(n) and create a matrix with size(k1)
287// colums and 2 rows, such that if [f1,f2] is a column of M then f1+f2 is one
288// of the polys of lin containing a pure degree 1 part and f1 is this part
289// interreduce this matrix (i.e. Gauss elimination on linear part, with rest
290// transformed accordingly).
291   int ii, kk;
292   if ( e!=0 )
293   {
294      ideal k1,k2;
295      int l = size(lin);
296      ideal lin1 = jet(lin,1) - jet(lin,0);  // part of pure degree 1
297      lin = lin - lin1;                      // rest, part of degree 1 substracted
298
299      for( ii=1; ii<=l; ii++ )
300      {
301         if( lin1[ii] != 0 )
302         {
303            kk = kk+1;
304            k1[kk] = lin1[ii];    // part of pure degree 1, renumbered
305            k2[kk] = lin[ii];     // rest of those polys which had a degree 1 part
306            lin[ii] = 0;
307         }
308      }
309
310      if( kk != 0 )
311      {
312         if( ordstr(P) != "c,dp(n)" )
313         {
314            execute("ring newP = ("+charstr(P)+"),("+varstr(P)+"),(c,dp);");
315            ideal k1 = imap(P,k1);
316            ideal k2 = imap(P,k2);
317            ideal lin = imap(P,lin);
318            ideal eva = imap(P,eva);
319            ideal sub = imap(P,sub);
320            ideal neva = imap(P,neva);
321         }
322         ideal k12 = k1,k2;
323         matrix M = matrix(k12,2,kk);
324        // M = interred(M);
325         l = ncols(M);
326         k1 = M[1,1..l];
327         k2 = M[2,1..l];
328         ideal kin = matrix(k1)+matrix(k2);
329         lin = simplify(lin,2);
330
331         l = size(kin);
332         poly p; map phi; ideal max;
333         for ( ii=1; ii<=n; ii++  )
334         {
335            for (kk=1; kk<=l; kk++ )
336            {
337               p = kin[kk]/var(ii);
338               if ( deg(p) == 0 )
339               {
340                  eva = eva+var(ii);
341                  neva[ii] = 0;
342                  sub = sub+kin[kk]/p;
343                  max = maxideal(1);
344                  max[ii] = var(ii) - (kin[kk]/p);
345                  phi = basering,max;
346                  lin = phi(lin);
347                  kin = simplify(phi(kin),2);
348                  l = size(kin);
349                 // ii=ii+1;
350                  break;
351               }
352            }
353         }
354         lin = kin+lin;
355      }
356   }
357   for( ii=1; ii<=size(lin); ii++ )
358   {
359      lin[ii] = cleardenom(lin[ii]);
360   }
361   if( defined(newP) )
362   {
363      setring P;
364      lin = imap(newP,lin);
365      eva = imap(newP,eva);
366      sub = imap(newP,sub);
367      neva = imap(newP,neva);
368   }
369   for( ii=1; ii<=n; ii++ )
370   {
371      for( kk=1; kk<=size(eva); kk++ )
372      {
373         if (phi[ii] == eva[kk] )
374         {  phi[ii] = eva[kk]-sub[kk]; break; }
375      }
376   }
377   map psi=P,phi;
378   ideal phi1=maxideal(1);
379   for(ii=1;ii<=size(eva);ii++){phi1=psi(phi1);}
380   L = lin, eva, sub, neva, phi1;
381  return(L);
382}
383example
384{ "EXAMPLE:"; echo = 2;
385   ring s=0,(x,y,z),dp;
386   ideal i =x2+y2,x2+y+1;
387   elimpart(i,3,0);
388   elimpart(i,3,1);
389}
390
391///////////////////////////////////////////////////////////////////////////////
392
393proc elimpartanyr (ideal i, list #)
394"USAGE:   elimpartanyr(i [,p,e] );  i=ideal, p=polynomial, e=integer@*
395         p: product of vars to be eliminated,@*
396         e =0: substitute from linear part of i (same as elimlinearpart)@*
397         e!=0: eliminate also by direct substitution@*
398         (default: p=product of all vars, e=1)
399RETURN:  list of 6 objects:
400  @format
401  [1]: (interreduced) ideal obtained by substituting from i those vars
402       appearing in p, which occur in the linear part of i (or which can
403       be expressed directly in the remaining variables, if e!=0)
404  [2]: ideal, variables which have been substituted
405  [3]: ideal, i-th element defines substitution of i-th var in [2]
406  [4]: ideal of variables of basering, substituted ones are set to 0
407  [5]: ideal, describing the map from the basering, say k[x(1..m)], to
408       itself onto k[..variables fom [4]..] and [1] is the image of i
409  [6]: int, # of vars considered for substitution (= # of factors of p)
410  @end format
411  The ideal i is generated by [1] and [3] in k[x(1..m)], the map [5]
412  maps [3] to 0, hence induces an isomorphism
413  @format
414            k[x(1..m)]/i -> k[..variables fom [4]..]/[1]
415  @end format
416NOTE:    the proc uses @code{execute} to create a ring with ordering dp and vars
417         placed correctly and then applies @code{elimpart}.
418EXAMPLE: example elimpartanyr; shows an example
419"
420{
421   def P = basering;
422   int j,n,e = 0,0,1;
423   poly p = product(maxideal(1));
424   if ( size(#)==1 ) { p=#[1]; }
425   if ( size(#)==2 ) { p=#[1]; e=#[2]; }
426   string a,b;
427   for ( j=1; j<=nvars(P); j++ )
428   {
429      if (deg(p/var(j))>=0) { a = a+varstr(j)+","; n = n+1; }
430      else { b = b+varstr(j)+","; }
431   }
432   if ( size(b) != 0 ) { b = b[1,size(b)-1]; }
433   else { a = a[1,size(a)-1]; }
434   execute("ring gnir ="+charstr(P)+",("+a+b+"),dp;");
435   ideal i = imap(P,i);
436   list L = elimpart(i,n,e)+list(n);
437   setring P;
438   list L = imap(gnir,L);
439   return(L);
440}
441example
442{ "EXAMPLE:"; echo = 2;
443   ring s=0,(x,y,z),dp;
444   ideal i = x3+y2+z,x2y2+z3,y+z+1;
445   elimpartanyr(i,z);
446}
447///////////////////////////////////////////////////////////////////////////////
448
449proc fastelim (ideal i, poly p, list #)
450"USAGE:   fastelim(i,p[h,o,a,b,e,m]); i=ideal, p=polynomial; h,o,a,b,e=integers
451          p: product of variables to be eliminated;@*
452  Optional parameters:
453  @format
454  - h !=0: use Hilbert-series driven std-basis computation
455  - o !=0: use proc @code{valvars} for a - hopefully - optimal ordering of vars
456  - a !=0: order vars to be eliminated w.r.t. increasing complexity
457  - b !=0: order vars not to be eliminated w.r.t. increasing complexity
458  - e !=0: use @code{elimpart} first to eliminate easy part
459  - m !=0: compute a minimal system of generators
460  @end format
461  (default: h,o,a,b,e,m = 0,1,0,0,0,0)
462RETURN:  ideal obtained from i by eliminating those variables, which occur in p
463EXAMPLE: example fastelim; shows an example.
464"
465{
466   def P = basering;
467   int h,o,a,b,e,m = 0,1,0,0,0,0;
468   if ( size(#) == 1 ) { h=#[1]; }
469   if ( size(#) == 2 ) { h=#[1]; o=#[2]; }
470   if ( size(#) == 3 ) { h=#[1]; o=#[2]; a=#[3]; }
471   if ( size(#) == 4 ) { h=#[1]; o=#[2]; a=#[3]; b=#[4];}
472   if ( size(#) == 5 ) { h=#[1]; o=#[2]; a=#[3]; b=#[4]; e=#[5]; }
473   if ( size(#) == 6 ) { h=#[1]; o=#[2]; a=#[3]; b=#[4]; e=#[5]; m=#[6]; }
474   list L = elimpartanyr(i,p,e);
475   poly q = product(L[2]);     //product of vars which are already eliminated
476   if ( q==0 ) { q=1; }
477   p = p/q;                    //product of vars which must still be eliminated
478   int nu = size(L[5])-size(L[2]);   //number of vars which must still be eliminated
479   if ( p==1 )                 //ready if no vars are left
480   {                           //compute minbase if 3-rd argument !=0
481      if ( m != 0 ) { L[1]=minbase(L[1]); }
482      return(L);
483   }
484//---------------- create new ring with remaining variables -------------------
485   string newvar = string(L[4]);
486   L = L[1],p;
487   execute("ring r1=("+charstr(P)+"),("+newvar+"),"+"dp;");
488   list L = imap(P,L);
489//------------------- find "best" ordering of variables  ----------------------
490   newvar = string(maxideal(1));
491   if ( o != 0 )
492   {
493      list ordevar = valvars(L[1],a,L[2],b);
494      intvec v = ordevar[1];
495      newvar=string(sort(maxideal(1),v)[1]);
496//------------ create new ring with "best" ordering of variables --------------
497      def r0=changevar(newvar);
498      setring r0;
499      list L = imap(r1,L);
500      kill r1;
501      def r1 = r0;
502      kill r0;
503   }
504//----------------- h==0: eliminate remaining vars directly -------------------
505   if ( h == 0 )
506   {
507      L[1] = eliminate(L[1],L[2]);
508      def r2 = r1;
509   }
510   else
511//------- h!=0: homogenize and compute Hilbert series using hilbvec ----------
512   {
513      intvec hi = hilbvec(L[1]);         // Hilbert series of i
514      execute("ring r2=("+charstr(P)+"),("+varstr(basering)+",@homo),dp;");
515      list L = imap(r1,L);
516      L[1] = homog(L[1],@homo);          // @homo = homogenizing var
517//---- use Hilbert-series to eliminate variables with Hilbert-driven std -----
518      L[1] = eliminate(L[1],L[2],hi);
519      L[1]=subst(L[1],@homo,1);          // dehomogenize by setting @homo=1
520   }
521   if ( m != 0 )                         // compute minbase
522   {
523      if ( #[1] != 0 ) { L[1] = minbase(L[1]); }
524   }
525   def id = L[1];
526   setring P;
527   return(imap(r2,id));
528}
529example
530{ "EXAMPLE:"; echo = 2;
531   ring s=31991,(e,f,x,y,z,t,u,v,w,a,b,c,d),dp;
532   ideal i = w2+f2-1, x2+t2+a2-1,  y2+u2+b2-1, z2+v2+c2-1,
533            d2+e2-1, f4+2u, wa+tf, xy+tu+ab;
534   fastelim(i,xytua,1,1);       //with hilb,valvars
535   fastelim(i,xytua,1,0,1);     //with hilb,minbase
536}
537///////////////////////////////////////////////////////////////////////////////
538
539proc faststd (@id, list #)
540"USAGE:   faststd(id [,\"hilb\",\"sort\",\"dec\",o,\"blocks\"]);
541         id=ideal/module, o=string (allowed:\"lp\",\"dp\",\"Dp\",\"ls\",
542         \"ds\",\"Ds\"),  \"hilb\",\"sort\",\"dec\",\"block\" options for
543         Hilbert-driven std, and the procedure sortandmap
544RETURN:  a ring R, in which an ideal STD_id is stored: @*
545         - the ring R differs from the active basering only in the choice
546         of monomial ordering and in the sorting of the variables.
547         - STD_id is a standard basis for the image (under imap) of the input
548         ideal/module id with respect to the new monomial ordering. @*
549NOTE:    Using the optional input parameters, we may modify the computations
550         performed: @*
551         - \"hilb\"  : use Hilbert-driven standard basis computation@*
552         - \"sort\"  : use 'sortandmap' for a best sorting of the variables@*
553         - \"dec\"   : order vars w.r.t. decreasing complexity (with \"sort\")@*
554         - \"block\" : create block ordering, each block having ordstr=o, s.t.
555                     vars of same complexity are in one block (with \"sort\")@*
556         - o       : defines the basic ordering of the resulting ring@*
557         [default: o=ordering of 1st block of basering (if allowed, else o=\"dp\"],
558                  \"sort\", if none of the optional parameters is given @*
559         This proc is only useful for hard problems where other methods fail.@*
560         \"hilb\" is useful for hard orderings (as \"lp\") or for characteristic 0,@*
561         it is correct for \"lp\",\"dp\",\"Dp\" (and for block orderings combining
562         these) but not for s-orderings or if the vars have different weights.@*
563         There seem to be only few cases in which \"dec\" is fast
564EXAMPLE: example faststd; shows an example.
565"
566{
567   def @P = basering;
568   int @h,@s,@n,@m,@ii = 0,0,0,0,0;
569   string @o,@va,@c = ordstr(basering),"","";
570//-------------------- prepare ordering and set options -----------------------
571   if ( @o[1]=="c" or @o[1]=="C")
572      {  @o = @o[3,2]; }
573   else
574      { @o = @o[1,2]; }
575   if( @o[1]!="d" and @o[1]!="D" and @o[1]!="l")
576      { @o="dp"; }
577
578   if (size(#) == 0 )
579      { @s = 1; }
580   for ( @ii=1; @ii<=size(#); @ii++ )
581   {
582      if ( typeof(#[@ii]) != "string" )
583      {
584         "// wrong syntax! type: help faststd";
585         return();
586      }
587      else
588      {
589         if ( #[@ii] == "hilb"  ) { @h = 1; }
590         if ( #[@ii] == "dec"   ) { @n = 1; }
591         if ( #[@ii] == "block" ) { @m = 1; }
592         if ( #[@ii] == "sort"  ) { @s = 1; }
593         if ( #[@ii]=="lp" or #[@ii]=="dp" or #[@ii]=="Dp" or #[@ii]=="ls"
594              or #[@ii]=="ds" or #[@ii]=="Ds" ) { @o = #[@ii]; }
595      }
596   }
597   if( voice==2 ) { "// chosen options, hilb sort dec block:",@h,@s,@n,@m; }
598
599//-------------------- nosort: create ring with new name ----------------------
600   if ( @s==0 )
601   {
602      execute("ring @S1 =("+charstr(@P)+"),("+varstr(@P)+"),("+@o+");");
603      def STD_id = imap(@P,@id);
604      if ( @h==0 ) { STD_id = std(STD_id); }
605   }
606
607//---------------------- no hilb: compute SB directly -------------------------
608   if ( @s != 0 and @h == 0 )
609   {
610      intvec getoption = option(get);
611      option(redSB);
612      @id = interred(sort(@id)[1]);
613      poly @p = product(maxideal(1),1..nvars(@P));
614      def @S1=sortandmap(@id,@n,@p,0,@o,@m);
615      setring @S1;
616      option(set,getoption);
617      def STD_id=imap(@S1,IMAG);
618      STD_id = std(STD_id);
619   }
620//------- hilb: homogenize and compute Hilbert-series using hilbvec -----------
621// this uses another standardbasis computation
622   if ( @h != 0 )
623   {
624      execute("ring @Q=("+charstr(@P)+"),("+varstr(@P)+",@homo),("+@o+");");
625      def @id = imap(@P,@id);
626      @id = homog(@id,@homo);               // @homo = homogenizing var
627      if ( @s != 0 )
628      {
629        intvec getoption = option(get);
630        option(redSB);
631        @id = interred(sort(@id)[1]);
632        poly @p = product(maxideal(1),1..(nvars(@Q)-1));
633        def @S1=sortandmap(@id,@n,@p,0,@o,@m);
634        setring @S1;
635        option(set,getoption);
636        kill @Q;
637        def @Q= basering;
638        def @id = IMAG;
639      }
640      intvec @hi;                     // encoding of Hilbert-series of i
641      @hi = hilbvec(@id);
642      //if ( @s!=0 ) { @hi = hilbvec(@id,"32003",ordstr(@Q)); }
643      //else { @hi = hilbvec(@id); }
644//-------------------------- use Hilbert-driven std --------------------------
645      @id = std(@id,@hi);
646      @id = subst(@id,@homo,1);             // dehomogenize by setting @homo=1
647      @va = varstr(@Q)[1,size(varstr(@Q))-6];
648      if ( @s!=0 )
649      {
650         @o = ordstr(@Q);
651         if ( @o[1]=="c" or @o[1]=="C") { @o = @o[1,size(@o)-6]; }
652         else { @o = @o[1,size(@o)-8] + @o[size(@o)-1,2]; }
653      }
654      kill @S1;
655      execute("ring @S1=("+charstr(@Q)+"),("+@va+"),("+@o+");");
656      def STD_id = imap(@Q,@id);
657   }
658   attrib(STD_id,"isSB",1);
659   export STD_id;
660   if (defined(IMAG)) { kill IMAG; }
661   setring @P;
662   dbprint(printlevel-voice+3,"
663// 'faststd' created a ring, in which an object STD_id is stored.
664// To access the object, type (if the name R was assigned to the return value):
665        setring R; STD_id; ");
666   return(@S1);
667}
668example
669{ "EXAMPLE:"; echo = 2;
670   system("--ticks-per-sec",100); // show time in 1/100 sec
671   ring s = 0,(e,f,x,y,z,t,u,v,w,a,b,c,d),(c,lp);
672   ideal i = w2+f2-1, x2+t2+a2-1,  y2+u2+b2-1, z2+v2+c2-1,
673            d2+e2-1, f4+2u, wa+tf, xy+tu+ab;
674   option(prot); timer=1;
675   int time = timer;
676   ideal j=std(i);
677   timer-time;
678   dim(j),mult(j);
679
680   time = timer;
681   def R=faststd(i);                      // use "best" ordering of vars
682   timer-time;
683   show(R);setring R;dim(STD_id),mult(STD_id);
684
685   setring s;kill R;time = timer;
686   def R=faststd(i,"hilb");                // hilb-std only
687   timer-time;
688   show(R);setring R;dim(STD_id),mult(STD_id);
689
690   setring s;kill R;time = timer;
691   def R=faststd(i,"hilb","sort");         // hilb-std,"best" ordering
692   timer-time;
693   show(R);setring R;dim(STD_id),mult(STD_id);
694
695   setring s;kill R;time = timer;
696   def R=faststd(i,"hilb","sort","block","dec"); // hilb-std,"best",blocks
697   timer-time;
698   show(R);setring R;dim(STD_id),mult(STD_id);
699
700   setring s;kill R;time = timer;
701   timer-time;time = timer;
702   def R=faststd(i,"sort","block","Dp"); //"best",decreasing,Dp-blocks
703   timer-time;
704   show(R);setring R;dim(STD_id),mult(STD_id);
705}
706///////////////////////////////////////////////////////////////////////////////
707
708proc findvars(id, list #)
709"USAGE:   findvars(id [,any] ); id=poly/ideal/vector/module/matrix, any=any type
710RETURN:  if no second argument is present: ideal of variables occuring in id,@*
711         if a second argument is given (of any type): list L with 4 entries:
712  @format
713  L[1]: ideal of variables occuring in id
714  L[2]: intvec of variables occuring in id
715  L[3]: ideal of variables not occuring in id
716  L[4]: intvec of variables not occuring in id
717  @end format
718EXAMPLE: example findvars; shows an example
719"
720{
721   int ii,n;
722   ideal found, notfound;
723   intvec f,nf;
724   n = nvars(basering);
725   ideal i = simplify(ideal(matrix(id)),10);
726   matrix M[ncols(i)][1] = i;
727   vector v = module(M)[1];
728   ideal max = maxideal(1);
729
730   for (ii=1; ii<=n; ii++)
731   {
732      if ( v != subst(v,var(ii),0) )
733      {
734         found = found+var(ii);
735         f = f,ii;
736      }
737      else
738      {
739         notfound = notfound+var(ii);
740         nf = nf,ii;
741      }
742   }
743   if ( size(f)>1 ) { f = f[2..size(f)]; }      //intvec of found vars
744   if ( size(nf)>1 ) { nf = nf[2..size(nf)]; }  //intvec of vars not found
745   if( size(#)==0 )  { return(found); }
746   if( size(#)!=0 )  { list L = found,f,notfound,nf; return(L); }
747}
748example
749{ "EXAMPLE:"; echo = 2;
750   ring s  = 0,(e,f,x,y,t,u,v,w,a,d),dp;
751   ideal i = w2+f2-1, x2+t2+a2-1;
752   findvars(i);
753   findvars(i,1);
754}
755///////////////////////////////////////////////////////////////////////////////
756
757proc hilbvec (@id, list #)
758"USAGE:   hilbvec(id[,c,o]); id=poly/ideal/vector/module/matrix, c,o=strings,@*
759          c=char, o=ordering used by @code{hilb}@*
760          (default: c=\"32003\", o=\"dp\")
761RETURN:  intvec of 1-st Hilbert-series of id, computed in char c and ordering o
762NOTE:    id must be homogeneous (i.e. all vars have weight 1)
763EXAMPLE: example hilbvec; shows an example
764"
765{
766   def @P = basering;
767   string @c,@o = "32003", "dp";
768   if ( size(#) == 1 ) {  @c = #[1]; }
769   if ( size(#) == 2 ) {  @c = #[1]; @o = #[2]; }
770   string @si = typeof(@id)+" @i = "+string(@id)+";";  //** weg
771   execute("ring @r=("+@c+"),("+varstr(basering)+"),("+@o+");");
772   //**def i = imap(P,@id);
773   execute(@si);                   //** weg
774   //show(basering);
775   @i = std(@i);
776   intvec @hi = hilb(@i,1);         // intvec of 1-st Hilbert-series of id
777   return(@hi);
778}
779example
780{ "EXAMPLE:"; echo = 2;
781   ring s   = 0,(e,f,x,y,z,t,u,v,w,a,b,c,d,H),dp;
782   ideal id = w2+f2-1, x2+t2+a2-1,  y2+u2+b2-1, z2+v2+c2-1,
783              d2+e2-1, f4+2u, wa+tf, xy+tu+ab;
784   id = homog(id,H);
785   hilbvec(id);
786}
787///////////////////////////////////////////////////////////////////////////////
788
789proc linearpart (id)
790"USAGE:   linearpart(id);  id=ideal/module
791RETURN:  generators of id of total degree <= 1
792EXAMPLE: example linearpart; shows an example
793"
794{
795   return(degreepart(id,0,1));
796}
797example
798{ "EXAMPLE:"; echo = 2;
799   ring r=0,(x,y,z),dp;
800   ideal i=1+x+x2+x3,3,x+3y+5z;
801   linearpart(i);
802   module m=[x,y,z],x*[x3,y2,z],[1,x2,z3,0,1];
803   show(linearpart(m));
804}
805///////////////////////////////////////////////////////////////////////////////
806
807proc tolessvars (id ,list #)
808"USAGE:   tolessvars(id [,s1,s2] ); id poly/ideal/vector/module/matrix,
809          s1=string (new ordering)@*
810          [default: s1=\"dp\" or \"ds\" depending on whether the first block
811          of the old ordering is a p- resp. an s-ordering)
812RETURN:  If id contains all vars of the basering: empty list. @*
813         Else: ring R with the same char as the basering, but possibly less
814         variables (only those variables which actually occur in id). In R
815         an object IMAG (image of id under imap) is stored.
816DISPLAY: If printlevel >=0, display ideal of vars, which have been ommitted
817         from the old ring.
818EXAMPLE: example tolessvars; shows an example
819"
820{
821//---------------- initialisation and check occurence of vars -----------------
822   int s,ii,n,fp,fs;
823   string s2,newvar;
824   int pr = printlevel-voice+3;  // p = printlevel+1 (default: p=1)
825   def P = basering;
826   s2 = ordstr(P);
827
828   list L = findvars(id,1);
829   newvar = string(L[1]);    // string of new variables
830   n = size(L[1]);           // number of new variables
831   if( n == 0 )
832   {
833      dbprint( pr,"","// no variable occured in "+typeof(id)+", no change of ring!");
834      return(id);
835   }
836   if( n == nvars(P) )
837   {
838     dbprint(printlevel-voice+3,"
839// All variables appear in input object.
840// empty list returned. ");
841     return(list());
842   }
843//----------------- prepare new ring, map to it and return --------------------
844   if ( size(#) == 0 )
845   {
846       fp = find(s2,"p");
847       fs = find(s2,"s");
848       if( fs==0 or (fs>=fp && fp!=0) ) { s2="dp"; }
849       else {  s2="ds"; }
850   }
851   if ( size(#) ==1 ) { s2=#[1]; }
852   dbprint( pr,"","// variables which did not occur:",L[3] );
853   execute("ring S1=("+charstr(P)+"),("+newvar+"),("+s2+");");
854   def IMAG = imap(P,id);
855   export IMAG;
856   dbprint(printlevel-voice+3,"
857// 'tolessvars' created a ring, in which an object IMAG is stored.
858// To access the object, type (if the name R was assigned to the return value):
859        setring R; IMAG; ");
860   return(S1);
861}
862example
863{ "EXAMPLE:"; echo = 2;
864   ring r  = 0,(x,y,z),dp;
865   ideal i = y2-x3,x-3,y-2x;
866   def R_r = tolessvars(i,"lp");
867   setring R_r;
868   show(basering);
869   IMAG;
870   kill R_r;
871}
872///////////////////////////////////////////////////////////////////////////////
873
874proc solvelinearpart (id,list #)
875"USAGE:   solvelinearpart(id [,n] );  id=ideal/module, n=integer,@*
876          (default: n=0)
877RETURN:  (interreduced) generators of id of degree <=1 in reduced triangular
878         form if n=0 [non-reduced triangular form if n!=0]
879ASSUME:  monomial ordering is a global ordering (p-ordering)
880NOTE:    may be used to solve a system of linear equations
881         see proc @code{gauss_row} from 'matrix.lib' for a different method
882WARNING: the result is very likely to be false for 'real' coefficients, use
884EXAMPLE: example solvelinearpart; shows an example
885"
886{
887   intvec getoption = option(get);
888   option(redSB);
889   if ( size(#)!=0 )
890   {
891      if(#[1]!=0) { option(noredSB); }
892   }
893   def lin = interred(degreepart(id,0,1));
894   if ( size(#)!=0 )
895   {
896      if(#[1]!=0)
897      {
898         return(lin);
899      }
900   }
901   option(set,getoption);
902   return(simplify(lin,1));
903}
904example
905{ "EXAMPLE:"; echo = 2;
906   // Solve the system of linear equations:
907   //         3x +   y +  z -  u = 2
908   //         3x +  8y + 6z - 7u = 1
909   //        14x + 10y + 6z - 7u = 0
910   //         7x +  4y + 3z - 3u = 3
911   ring r = 0,(x,y,z,u),lp;
912   ideal i= 3x +   y +  z -  u,
913           13x +  8y + 6z - 7u,
914           14x + 10y + 6z - 7u,
915            7x +  4y + 3z - 3u;
916   ideal j= 2,1,0,3;
917   j = i-j;                        // difference of 1x4 matrices
918                                   // compute reduced triangular form, setting
919   solvelinearpart(j);             // the RHS equal 0 gives the solutions!
920   solvelinearpart(j,1); "";       // triangular form, not reduced
921}
922///////////////////////////////////////////////////////////////////////////////
923
924proc sortandmap (@id, list #)
925"USAGE:   sortandmap(id [,n1,p1,n2,p2...,o1,m1,o2,m2...]);@*
926         id=poly/ideal/vector/module,@*
927         p1,p2,...= polynomials (product of variables),@*
928         n1,n2,...= integers,@*
929         o1,o2,...= strings,@*
930         m1,m2,...= integers@*
931         (default: p1=product of all vars, n1=0, o1=\"dp\",m1=0)
932         the last pi (containing the remaining vars) may be omitted
933RETURN:  a ring R, in which a poly/ideal/vector/module IMAG is stored: @*
934         - the ring R differs from the active basering only in the choice
935         of monomial ordering and in the sorting of the variables.@*
936         - IMAG is the image (under imap) of the input ideal/module id @*
937         The new monomial ordering and sorting of vars is as follows:
938  @format
939  - each block of vars occuring in pi is sorted w.r.t. its complexity in id,
940  - ni controls the sorting in i-th block (= vars occuring in pi):
941    ni=0 (resp.!=0) means that least complex (resp. most complex) vars come
942    first @*
943  - oi and mi define the monomial ordering of the i-th block:
944    if mi =0, oi=ordstr(i-th block)
945    if mi!=0, the ordering of the i-th block itself is a blockordering,
946      each subblock having ordstr=oi, such that vars of same complexity are
947      in one block
948  @end format
949         Note that only simple ordstrings oi are allowed: \"lp\",\"dp\",\"Dp\",
950         \"ls\",\"ds\",\"Ds\". @*
951NOTE:    We define a variable x to be more complex than y (with respect to id)
952         if val(x) > val(y) lexicographically, where val(x) denotes the
953         valuation vector of x:@*
954         consider id as list of polynomials in x with coefficients in the
955         remaining variables. Then:@*
956         val(x) = (maximal occuring power of x,  # of all monomials in leading
957         coefficient, # of all monomials in coefficient of next smaller power
958         of x,...).
959EXAMPLE: example sortandmap; shows an example
960"
961{
962   def @P = basering;
963   int @ii,@jj;
964   intvec @v;
965   string @o;
966//----------------- find o in # and split # into 2 lists ---------------------
967   # = # +list("dp",0);
968   for ( @ii=1; @ii<=size(#); @ii++)
969   {
970      if ( typeof(#[@ii])=="string" )  break;
971   }
972   if ( @ii==1 ) { list @L1 = list(); }
973   else { list @L1 = #[1..@ii-1]; }
974   list @L2 = #[@ii..size(#)];
975   list @L = sortvars(@id,@L1);
976   string @va = string(@L[1]);
977   list @l = @L[2];   //e.g. @l[4]=intvec describing permutation of 1-st block
978//----------------- construct correct ordering with oi and mi ----------------
979   for ( @ii=4; @ii<=size(@l); @ii=@ii+4 )
980   {
981      @L2=@L2+list("dp",0);
982      if ( @L2[@ii/2] != 0)
983      {
984         @v = @l[@ii];
985         for ( @jj=1; @jj<=size(@v); @jj++ )
986         {
987           @o = @o+@L2[@ii/2 -1]+"("+string(@v[@jj])+"),";
988         }
989      }
990      else
991      {
992         @o = @o+@L2[@ii/2 -1]+"("+string(size(@l[@ii/2]))+"),";
993      }
994   }
995   @o=@o[1..size(@o)-1];
996   execute("ring @S1 =("+charstr(@P)+"),("+@va+"),("+@o+");");
997   def IMAG = imap(@P,@id);
998   export IMAG;
999   dbprint(printlevel-voice+3,"
1000// 'sortandmap' created a ring, in which an object IMAG is stored.
1001// To access the object, type (if the name R was assigned to the return value):
1002        setring R; IMAG; ");
1003   return(@S1);
1004}
1005example
1006{ "EXAMPLE:"; echo = 2;
1007   ring s = 32003,(x,y,z),dp;
1008   ideal i=x3+y2,xz+z2;
1009   def R_r=sortandmap(i);
1010   show(R_r);
1011   setring R_r; IMAG;
1012   kill R_r; setring s;
1013   def R_r=sortandmap(i,1,xy,0,z,0,"ds",0,"lp",0);
1014   show(R_r);
1015   setring R_r; IMAG;
1016   kill R_r;
1017}
1018///////////////////////////////////////////////////////////////////////////////
1019
1020proc sortvars (id, list #)
1021"USAGE:   sortvars(id[,n1,p1,n2,p2,...]);@*
1022         id=poly/ideal/vector/module,@*
1023         p1,p2,...= polynomials (product of vars),@*
1024         n1,n2,...=integers@*
1025         (default: p1=product of all vars, n1=0)
1026         the last pi (containing the remaining vars) may be omitted
1027COMPUTE: sort variables with respect to their complexity in id
1028RETURN:  list of two elements, an ideal and a list:
1029  @format
1030  [1]: ideal, variables of basering sorted w.r.t their complexity in id
1031       ni controls the ordering in i-th block (= vars occuring in pi):
1032       ni=0 (resp.!=0) means that less (resp. more) complex vars come first
1033  [2]: a list with 4 entries for each pi:
1034       ideal ai : vars of pi in correct order,
1035       intvec vi: permutation vector describing the ordering in ai,
1036       intmat Mi: valuation matrix of ai, the columns of Mi being the
1037                  valuation vectors of the vars in ai
1038       intvec wi: size of 1-st, 2-nd,... block of identical columns of Mi
1039                  (vars with same valuation)
1040  @end format
1041NOTE:    We define a variable x to be more complex than y (with respect to id)
1042         if val(x) > val(y) lexicographically, where val(x) denotes the
1043         valuation vector of x:@*
1044         consider id as list of polynomials in x with coefficients in the
1045         remaining variables. Then:@*
1046         val(x) = (maximal occuring power of x,  # of all monomials in leading
1047         coefficient, # of all monomials in coefficient of next smaller power
1048         of x,...).
1049EXAMPLE: example sortvars; shows an example
1050"
1051{
1052   int ii,jj,n,s;
1053   list L = valvars(id,#);
1054   list L2, L3 = L[2], L[3];
1055   list K; intmat M; intvec v1,v2,w;
1056   ideal i = sort(maxideal(1),L[1])[1];
1057   for ( ii=1; ii<=size(L2); ii++ )
1058   {
1059      M = transpose(L3[2*ii]);
1060      M = M[L2[ii],1..nrows(L3[2*ii])];
1061      w = 0; s = 0;
1062      for ( jj=1; jj<=nrows(M)-1; jj++ )
1063      {
1064         v1 = M[jj,1..ncols(M)];
1065         v2 = M[jj+1,1..ncols(M)];
1066         if ( v1 != v2 ) { n=jj-s; s=s+n; w = w,n; }
1067      }
1068      w=w,nrows(M)-s; w=w[2..size(w)];
1069      K = K+sort(L3[2*ii-1],L2[ii])+list(transpose(M))+list(w);
1070   }
1071   L = i,K;
1072   return(L);
1073}
1074example
1075{ "EXAMPLE:"; echo = 2;
1076   ring s=0,(x,y,z,w),dp;
1077   ideal i = x3+y2+yw2,xz+z2,xyz-w2;
1078   sortvars(i,0,xy,1,zw);
1079}
1080///////////////////////////////////////////////////////////////////////////////
1081
1082proc valvars (id, list #)
1083"USAGE:   valvars(id[,n1,p1,n2,p2,...]);@*
1084         id=poly/ideal/vector/module,@*
1085         p1,p2,...= polynomials (product of vars),@*
1086         n1,n2,...= integers,
1087
1088         ni controls the ordering of vars occuring in pi: ni=0 (resp.!=0) means
1089         that less (resp. more) complex vars come first@*
1090         (default: p1=product of all vars, n1=0)
1091         the last pi (containing the remaining vars) may be omitted
1092COMPUTE: valuation (complexity) of variables with respect to id.@*
1093         ni controls the ordering of vars occuring in pi:@*
1094         ni=0 (resp.!=0) means that less (resp. more) complex vars come first.
1095RETURN:  list with 3 entries:
1096  @format
1097  [1]: intvec, say v, describing the permutation such that the permuted
1098       ringvariables are ordered with respect to their complexity in id
1099  [2]: list of intvecs, i-th intvec, say v(i) describing permutation
1100       of vars in a(i) such that v=v(1),v(2),...
1101  [3]: list of ideals and intmat's, say a(i) and M(i), where
1102       a(i): factors of pi,
1103       M(i): valuation matrix of a(i), such that the j-th column of M(i)
1104             is the valuation vector of j-th generator of a(i)
1105         @end format
1106NOTE:    Use @code{sortvars} in order to actually sort the variables!
1107         We define a variable x to be more complex than y (with respect to id)
1108         if val(x) > val(y) lexicographically, where val(x) denotes the
1109         valuation vector of x:@*
1110         consider id as list of polynomials in x with coefficients in the
1111         remaining variables. Then:@*
1112         val(x) = (maximal occuring power of x,  # of all monomials in leading
1113         coefficient, # of all monomials in coefficient of next smaller power
1114         of x,...).
1115EXAMPLE: example valvars; shows an example
1116"
1117{
1118//---------------------------- initialization ---------------------------------
1119   int ii,jj,kk,n;
1120   list L;                    // list of valuation vectors in one block
1121   intvec vec;                // describes permutation of vars (in one block)
1122   list blockvec;             // i-th element = vec of i-th block
1123   intvec varvec;             // result intvector
1124   list Li;                   // result list of ideals
1125   list LM;                   // result list of intmat's
1126   intvec v,w,s;              // w valuation vector for one variable
1127   matrix C;                  // coefficient matrix for different variables
1128   ideal i = simplify(ideal(matrix(id)),10);
1129
1130//---- for each pii in # create ideal a(ii) intvec v(ii) and list L(ii) -------
1131// a(ii) = ideal of vars in product, v(ii)[j]=k <=> a(ii)[j]=var(k)
1132
1133   v = 1..nvars(basering);
1134   int l = size(#);
1135   if ( l >= 2 )
1136   {
1137      ideal m=maxideal(1);
1138      for ( ii=2; ii<=l; ii=ii+2 )
1139      {
1140         int n(ii) = #[ii-1];
1141         ideal a(ii);
1142         intvec v(ii);
1143         for ( jj=1; jj<=nvars(basering); jj++ )
1144         {
1145            if ( #[ii]/var(jj) != 0)
1146            {
1147               a(ii) = a(ii) + var(jj);
1148               v(ii)=v(ii),jj;
1149               m[jj]=0;
1150               v[jj]=0;
1151            }
1152         }
1153         v(ii)=v(ii)[2..size(v(ii))];
1154      }
1155      if ( size(m)!=0 )
1156      {
1157         l = 2*(l/2)+2;
1158         ideal a(l) = simplify(m,2);
1159         intvec v(l) = compress(v);
1160         int n(l);
1161         if ( size(#)==l-1 ) { n(l) = #[l-1]; }
1162      }
1163   }
1164   else
1165   {
1166      l = 2;
1167      ideal a(2) = maxideal(1);
1168      intvec v(2) = v;
1169      int n(2);
1170      if ( size(#)==1 ) { n(2) = #[1]; }
1171   }
1172//------------- start loop to order variables in each a(ii) -------------------
1173
1174   for ( kk=2; kk<=l; kk=kk+2 )
1175   {
1176      L = list();
1177      n = 0;
1178//---------------- get valuation of all variables in a(kk) --------------------
1179      for ( ii=1; ii<=size(a(kk)); ii++ )
1180      {
1181         C = coeffs(i,a(kk)[ii]);
1182         w = nrows(C); // =(maximal occuring power of a(kk)[ii])+1
1183         for ( jj=w[1]; jj>1; jj-- )
1184         {
1185            s = size(C[jj,1..ncols(C)]);
1186            w[w[1]-jj+2] = sum(s);
1187         }
1188         // w[1] should represent the maximal occuring power of a(kk)[ii] so it
1189         // has to be decreased by 1 since otherwise the constant term is also
1190         // counted
1191         w[1]=w[1]-1;
1192
1193         L[ii]=w;
1194         n = size(w)*(size(w) > n) + n*(size(w) <= n);
1195      }
1196      intmat M(kk)[size(a(kk))][n];
1197      for ( ii=1; ii<=size(a(kk)); ii++ )
1198      {
1199         if ( n==1 ) { w = L[ii]; M(kk)[ii,1] = w[1]; }
1200         else  { M(kk)[ii,1..n] = L[ii]; }
1201      }
1202      LM[kk-1] = a(kk);
1203      LM[kk] = transpose(compress(M(kk)));
1204//------------------- compare valuation and insert in vec ---------------------
1205      vec = sort(L)[2];
1206      if ( n(kk) != 0 ) { vec = vec[size(vec)..1]; }
1207      blockvec[kk/2] = vec;
1208      vec = sort(v(kk),vec)[1];
1209      varvec = varvec,vec;
1210   }
1211   varvec = varvec[2..size(varvec)];
1212   list result = varvec,blockvec,LM;
1213   return(result);
1214}
1215example
1216{ "EXAMPLE:"; echo = 2;
1217   ring s=0,(x,y,z,a,b),dp;
1218   ideal i=ax2+ay3-b2x,abz+by2;
1219   valvars (i,0,xyz);
1220}
1221///////////////////////////////////////////////////////////////////////////////
1222proc idealSplit(ideal I,list #)
1223"USAGE:  idealSplit(id,timeF,timeS);  id ideal and optional
1224         timeF ,timeS integers to bound the time which can be used
1225         for factorization resp. standard basis computation
1226RETURN:  a list of ideals such that their intersection
1227         has the same radical as id
1228EXAMPLE: example idealSplit; shows an example
1229"
1230{
1231   option(redSB);
1232   int j,k,e;
1233   int i=1;
1234   int l=attrib(I,"isSB");
1235   ideal J;
1236   int timeF;
1237   int timeS;
1238   list re,fac,te;
1239
1240   if(size(#)==1)
1241   {
1242     if(typeof(#[1])=="ideal")
1243     {
1244        re=#;
1245     }
1246     else
1247     {
1248       timeF=#[1];
1249     }
1250   }
1251   if(size(#)==2)
1252   {
1253     if(typeof(#[1])=="list")
1254     {
1255        re=#[1];
1256        timeF=#[2];
1257     }
1258     else
1259     {
1260       timeF=#[1];
1261       timeS=#[2];
1262     }
1263   }
1264   if(size(#)==3){re=#[1];timeF=#[2];timeS=#[3];}
1265
1266   fac=timeFactorize(I[1],timeF);
1267
1268   while((size(fac[1])==2)&&(i<size(I)))
1269   {
1270      i++;
1271      fac=timeFactorize(I[i],timeF);
1272   }
1273   if(size(fac[1])>2)
1274   {
1275      for(j=2;j<=size(fac[1]);j++)
1276      {
1277         I[i]=fac[1][j];
1278         attrib(I,"isSB",1);
1279         e=1;
1280         k=0;
1281         while(k<size(re))
1282         {
1283            k++;
1284            if(size(reduce(re[k],I))==0){e=0;break;}
1285            attrib(re[k],"isSB",1);
1286            if(size(reduce(I,re[k]))==0){re=delete(re,k);k--;}
1287         }
1288         if(e)
1289         {
1290            if(l)
1291            {
1292               J=I;
1293               J[i]=0;
1294               J=simplify(J,2);
1295               attrib(J,"isSB",1);
1296               re=idealSplit(std(J,fac[1][j]),re,timeF,timeS);
1297            }
1298            else
1299            {
1300               re=idealSplit(timeStd(I,timeS),re,timeF,timeS);
1301            }
1302         }
1303      }
1304      return(re);
1305   }
1306   J=timeStd(I,timeS);
1307   attrib(I,"isSB",1);
1308   if(size(reduce(J,I))==0){return(re+list(I));}
1309   return(re+idealSplit(J,re,timeF,timeS));
1310}
1311example
1312{ "EXAMPLE:"; echo = 2;
1313   ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
1314   ideal i=
1315   bv+su,
1316   bw+tu,
1317   sw+tv,
1318   by+sx,
1319   bz+tx,
1320   sz+ty,
1321   uy+vx,
1322   uz+wx,
1323   vz+wy,
1324   bvz;
1325   idealSplit(i);
1326}
1327///////////////////////////////////////////////////////////////////////////////
1328proc idealSimplify(ideal J,list #)
1329"USAGE:  idealSimplify(id);  id ideal
1330RETURN:  ideal I = eliminate(Id,m) m is a product of variables
1331         which are only linearly involved in the generators of id
1332EXAMPLE: example idealSimplify; shows an example
1333"
1334{
1335   ideal I=J;
1336   if(size(#)!=0){I=#[1];}
1337   def R=basering;
1338   matrix M=jacob(I);
1339   ideal ma=maxideal(1);
1340   int i,j,k;
1341   map phi;
1342
1343   for(i=1;i<=nrows(M);i++)
1344   {
1345      for(j=1;j<=ncols(M);j++)
1346      {
1347         if(deg(M[i,j])==0)
1348         {
1349            ma[j]=(-1/M[i,j])*(I[i]-M[i,j]*var(j));
1350            phi=R,ma;
1351            I=phi(I);
1352            J=phi(J);
1353            for(k=1;k<=ncols(I);k++){I[k]=cleardenom(I[k]);}
1354            M=jacob(I);
1355         }
1356      }
1357   }
1358   J=simplify(J,2);
1359   for(i=1;i<=size(J);i++){J[i]=cleardenom(J[i]);}
1360   return(J);
1361}
1362example
1363{ "EXAMPLE:"; echo = 2;
1364   ring r=0,(x,y,z,w,t),dp;
1365   ideal i=
1366   t,
1367   x3+y2+2z,
1368   x2+3y,
1369   x2+y2+z2,
1370   w2+z;
1371   ideal j=idealSimplify(i);
1372   ideal k=eliminate(i,zyt);
1373   reduce(k,std(j));
1374   reduce(j,std(k));
1375}
1376
1377///////////////////////////////////////////////////////////////////////////////
1378
1379/*
1380
1381 ring s=31991,(e,f,x,y,z,t,u,v,w,a,b,c,d),dp;
1382 ring s=31991,(x,y,z,t,u,v,w,a,b,c,d,f,e,h),dp; //standard
1383 ring s1=31991,(y,u,b,c,a,z,t,x,v,d,w,e,f,h),dp; //gut
1384v;
138513,12,11,10,8,7,6,5,4,3,2,1,9,14
1386print(matrix(sort(maxideal(1),v)));
1387f,e,w,d,x,t,z,a,c,b,u,y,v,h
1388print(matrix(maxideal(1)));
1389y,u,b,c,a,z,t,x,v,d,w,e,f,h
1390v0;
139114,9,12,11,10,8,7,6,5,4,3,2,1,13
1392print(matrix(sort(maxideal(1),v0)));
1393h,v,e,w,d,x,t,z,a,c,b,u,y,f
1394v1;v2;
13959,12,11,10,8,7,6,5,4,3,2,1,13,14
139613,12,11,10,8,7,6,5,4,3,2,1,9,14
1397
1398Ev. Gute Ordnung fuer i:
1399========================
1401mit ar=(ar1,...,ark), k=size(i)
1402    arj in K[..x^..]
1403d=deg_x(i) := max{deg_x(i[k]) | k=1..size(i)}
1405x>y  <==
1406  1. deg_x(i)>deg_y(i)
1407  2. "=" in 1. und size_x lexikographisch
1408
1409hier im Beispiel:
1410f: 5,1,0,1,2
1411
1412u: 3,1,4
1413
1414y: 3,1,3
1415b: 3,1,3
1416c: 3,1,3
1417a: 3,1,3
1418z: 3,1,3
1419t: 3,1,3
1420
1421x: 3,1,2
1422v: 3,1,2
1423d: 3,1,2
1424w: 3,1,2
1425e: 3,1,2
1426probier mal:
1427 ring s=31991,(f,u,y,z,t,a,b,c,v,w,d,e,h),dp; //standard
1428
1429*/
Note: See TracBrowser for help on using the repository browser.