source: git/Singular/LIB/presolve.lib @ 3ea5cee

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