source: git/Singular/LIB/presolve.lib @ 50cbdc

//last change: 13.02.2001 (Eric Westenberger)
///////////////////////////////////////////////////////////////////////////////
version="$Id: presolve.lib,v 1.18 2001-08-27 14:47:57 Singular Exp $";
category="Symbolic-numerical solving";
info="
LIBRARY:  presolve.lib     Pre-Solving of Polynomial Equations
AUTHOR:   Gert-Martin Greuel, email: greuel@mathematik.uni-kl.de,

PROCEDURES:
 degreepart(id,d1,d2);  elements of id of total degree >= d1 and <= d2
 elimlinearpart(id);    linear part eliminated from id
 elimpart(id[,n]);      partial elimination of vars [among first n vars]
 elimpartanyr(i,p);     factors of p partially eliminated from i in any ring
 fastelim(i,p[..]);     fast elimination of factors of p from i [options]
 findvars(id[..]);      ideal of variables occuring in id [more information]
 hilbvec(id[,c,o]);     intvec of Hilberseries of id [in char c and ord o]
 linearpart(id);        elements of id of total degree <=1
 tolessvars(id[,]);     maps id to new basering having only vars occuring in id
 solvelinearpart(id);   reduced std-basis of linear part of id
 sortandmap(id,s1,s2);  map to new basering with vars sorted w.r.t. complexity
 sortvars(id[n1,p1..]); sort vars w.r.t. complexity in id [different blocks]
 valvars(id[..]);       valuation of vars w.r.t. to their complexity in id
           (parameters in square brackets [] are optional)
";

LIB "inout.lib";
LIB "general.lib";
LIB "matrix.lib";
LIB "ring.lib";
LIB "elim.lib";
///////////////////////////////////////////////////////////////////////////////

proc degreepart (id,int d1,int d2,list #)
"USAGE:   degreepart(id,d1,d2[,v]);  id=ideal/module, d1,d1=integers, v=intvec
RETURN:  generators of id of [v-weighted] total degree >= d1 and <= d2
         (default: v = 1,...,1)
EXAMPLE: example degreepart; shows an example
"
{
   def dpart = id;
   int s,ii = size(id),0;
   if ( size(#)==0 )
   {
      for ( ii=1; ii<=s; ii=ii+1 )
      {
         dpart[ii] = (jet(id[ii],d1-1)==0)*(id[ii]==jet(id[ii],d2))*id[ii];
      }
   }
   else
   {
      for ( ii=1; ii<=s; ii=ii+1 )
      {
         dpart[ii]=(jet(id[ii],d1-1,#[1])==0)*(id[ii]==jet(id[ii],d2,#[1]))*id[ii];
      }
   }
   return(simplify(dpart,2));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   ideal i=1+x+x2+x3+x4,3,xz+y3+z8;
   degreepart(i,0,4);
   module m=[x,y,z],x*[x3,y2,z],[1,x2,z3,0,1];
   intvec v=2,3,6;
   show(degreepart(m,8,8,v));
}
///////////////////////////////////////////////////////////////////////////////

proc elimlinearpart (ideal i,list #)
"USAGE:   elimlinearpart(i[,n]);  i=ideal, n=integer,@*
          default: n=nvars(basering)
RETURN:   list L with 5 entries:
  @format
  L[1]: (interreduced) ideal obtained from i by substituing
        from the first n variables those, which appear in a linear part
        of i, by putting this part into triangular form
  L[2]: ideal of variables which have been substituted
  L[3]: ideal, j-th element defines substitution of j-th var in [2]
  L[4]: ideal of variables of basering, eliminated ones are set to 0
  L[5]: ideal, describing the map from the basering to itself such that
        L[1] is the image of i
  @end format
NOTE:    the procedure does always interreduce the ideal i internally w.r.t.
         ordering dp.
EXAMPLE: example elimlinearpart; shows an example
"
{
   int ii,n,fi,k;
   string o, newo;
   intvec getoption = option(get);
   option(redSB);
   def P = basering;
   n = nvars(P);
//--------------- replace ordering by dp-ordering if necessary ----------------
   o = "dp("+string(n)+")";
   fi = find(ordstr(P),o);
   if( fi == 0 or find(ordstr(P),"a") != 0 )
   {
      execute("ring newP = ("+charstr(P)+"),("+varstr(P)+"),dp;");
      ideal i = imap(P,i);
   }
//---------------------------------- start ------------------------------------
   if ( size(#)!=0 ) {  n=#[1]; }
   ideal maxi,rest = maxideal(1),0;
   if ( n < nvars(P) ) { rest = maxi[n+1..nvars(P)]; }
   attrib(rest,"isSB",1);
//-------------------- interreduce and find linear part  ----------------------
// interred does the only real work. Because of ordering dp the linear part is
// within the first elements after interreduction
// **: perhaps Bareiss to constant matrix of linear part instead of interred,
// and/or for big systems, interred only those generators of id
// which do not contain elements not to be eliminated

   ideal id = interred(i);

   if(id[1] == 1 )                          //special case of unit ideal
   {
     setring P;
     list L = ideal(1), ideal(0), ideal(0), maxideal(1), maxideal(1);
     option(set,getoption);
     return(L);
   }

   for ( ii=1; ii<=size(id); ii++ )
   {
      if( deg(id[ii]) > 1) { break; }
      k=ii;
   }
   if( k == 0 )       { ideal lin; }
   else               { ideal lin = id[1..k];}
   if( k < size(id) ) { id = id[k+1..size(id)]; }
   else               { id = 0; }
//----- remove generators from lin containing vars not to be eliminated  ------
   if ( n < nvars(P) )
   {
      for ( ii=1; ii<=size(lin); ii++ )
      {
         if ( reduce(lead(lin[ii]),rest) == 0 )
         {
            id=lin[ii],id;
            lin[ii] = 0;
         }
      }
   }
   lin = simplify(lin,2);
   attrib(lin,"isSB",1);
   ideal eva = lead(lin);
   attrib(eva,"isSB",1);
   ideal neva = reduce(maxideal(1),eva);
//------------------ go back to original ring end return  ---------------------
   if( fi == 0 or find(ordstr(P),"a") != 0 )
   {
      setring P;
      ideal id = imap(newP,id);
      ideal eva = imap(newP,eva);
      ideal lin = imap(newP,lin);
      ideal neva = imap(newP,neva);
   }

   eva = eva[ncols(eva)..1];  // sorting according to variables in basering
   lin = lin[ncols(lin)..1];
   ideal phi= neva;
   k = 1;
   for( ii=1; ii<=n; ii++ )
   {
      if( neva[ii] == 0 )
      {
         phi[ii] = eva[k]-lin[k];
         k=k+1;
      }
   }
   list L = id, eva, lin, neva, phi;
   option(set,getoption);
   return(L);
}
example
{ "EXAMPLE:"; echo = 2;
   ring s=0,(x,y,z),dp;
   ideal i = x3+y2+z,x2y2+z3,y+z+1;
  elimlinearpart(i);
}
///////////////////////////////////////////////////////////////////////////////

proc elimpart (ideal i,list #)
"USAGE:   elimpart(i [,n,e] );  i=ideal, n,e=integers
         n   : only the first n vars are considered for substitution,@*
         e =0: substitute from linear part of i (same as elimlinearpart)@*
         e!=0: eliminate also by direct substitution@*
         (default: n = nvars(basering), e = 1)
RETURN:  list of 5 objects:
  @format
  [1]: ideal obtained by substituting from the first n variables those
       from i, which appear in the linear part of i (or, if e!=0, which
       can be expressed directly in the remaining vars)
  [2]: ideal, variables which have been substituted
  [3]: ideal, i-th element defines substitution of i-th var in [2]
  [4]: ideal of variables of basering, substituted ones are set to 0
  [5]: ideal, describing the map from the basering, say k[x(1..m)], to
       itself onto k[..variables fom [4]..] and [1] is the image of i
  @end format
  The ideal i is generated by [1] and [3] in k[x(1..m)], the map [5]
  maps [3] to 0, hence induces an isomorphism
  @format
            k[x(1..m)]/i -> k[..variables fom [4]..]/[1]
  @end format
NOTE:    If the basering has ordering (c,dp), this is faster for big ideals,
         since it avoids internal ring change and mapping.
EXAMPLE: example elimpart; shows an example
"
{
   def P = basering;
   int n,e = nvars(P),1;
   if ( size(#)==1 ) {  n=#[1]; }
   if ( size(#)==2 ) {  n=#[1]; e=#[2];}
//----------- interreduce linear part with proc elimlinearpart ----------------
// lin = ideal after interreduction of linear part
// eva = eliminated (substituted) variables
// sub = polynomials defining substitution
// neva= not eliminated variables
// phi = map describing substitution
   list L = elimlinearpart(i,n);
   ideal lin, eva, sub, neva, phi = L[1], L[2], L[3], L[4], L[5];
//-------- direct substitution of variables if possible and if e!=0 -----------
// first find terms lin1 of pure degree 1 in each poly of lin
// k1 = pure degree 1 part, k1+k2 = those polys of lin which contained a pure
// degree 1 part.
// Then go to ring newP with ordering c,dp(n) and create a matrix with size(k1)
// colums and 2 rows, such that if [f1,f2] is a column of M then f1+f2 is one
// of the polys of lin containing a pure degree 1 part and f1 is this part
// interreduce this matrix (i.e. Gauss elimination on linear part, with rest
// transformed accordingly).
   int ii, kk;
   if ( e!=0 )
   {
      ideal k1,k2;
      int l = size(lin);
      ideal lin1 = jet(lin,1) - jet(lin,0);  // part of pure degree 1
      lin = lin - lin1;                      // rest, part of degree 1 substracted

      for( ii=1; ii<=l; ii++ )
      {
         if( lin1[ii] != 0 )
         {
            kk = kk+1;
            k1[kk] = lin1[ii];    // part of pure degree 1, renumbered
            k2[kk] = lin[ii];     // rest of those polys which had a degree 1 part
            lin[ii] = 0;
         }
      }

      if( kk != 0 )
      {
         if( ordstr(P) != "c,dp(n)" )
         {
            execute("ring newP = ("+charstr(P)+"),("+varstr(P)+"),(c,dp);");
            ideal k1 = imap(P,k1);
            ideal k2 = imap(P,k2);
            ideal lin = imap(P,lin);
            ideal eva = imap(P,eva);
            ideal sub = imap(P,sub);
            ideal neva = imap(P,neva);
         }
         ideal k12 = k1,k2;
         matrix M = matrix(k12,2,kk);
        // M = interred(M);
         l = ncols(M);
         k1 = M[1,1..l];
         k2 = M[2,1..l];
         ideal kin = matrix(k1)+matrix(k2);
         lin = simplify(lin,2);

         l = size(kin);
         poly p; map phi; ideal max;
         for ( ii=1; ii<=n; ii++  )
         {
            for (kk=1; kk<=l; kk++ )
            {
               p = kin[kk]/var(ii);
               if ( deg(p) == 0 )
               {
                  eva = eva+var(ii);
                  neva[ii] = 0;
                  sub = sub+kin[kk]/p;
                  max = maxideal(1);
                  max[ii] = var(ii) - (kin[kk]/p);
                  phi = basering,max;
                  lin = phi(lin);
                  kin = simplify(phi(kin),2);
                  l = size(kin);
                  ii=ii+1;
                  break;
               }
            }
         }
         lin = kin+lin;
      }
   }
   for( ii=1; ii<=size(lin); ii++ )
   {
      lin[ii] = cleardenom(lin[ii]);
   }
   if( defined(newP) )
   {
      setring P;
      lin = imap(newP,lin);
      eva = imap(newP,eva);
      sub = imap(newP,sub);
      neva = imap(newP,neva);
   }
   for( ii=1; ii<=n; ii++ )
   {
      for( kk=1; kk<=size(eva); kk++ )
      {
         if (phi[ii] == eva[kk] )
         {  phi[ii] = eva[kk]-sub[kk]; break; }
      }
   }
   L = lin, eva, sub, neva, phi;
  return(L);
}
example
{ "EXAMPLE:"; echo = 2;
   ring s=0,(x,y,z),dp;
   ideal i =x2+y2,x2+y+1;
   elimpart(i,3,0);
   elimpart(i,3,1);
}

///////////////////////////////////////////////////////////////////////////////

proc elimpartanyr (ideal i, list #)
"USAGE:   elimpartanyr(i [,p,e] );  i=ideal, p=polynomial, e=integer@*
         p: product of vars to be eliminated,@*
         e =0: substitute from linear part of i (same as elimlinearpart)@*
         e!=0: eliminate also by direct substitution@*
         (default: p=product of all vars, e=1)
RETURN:  list of 6 objects:
  @format
  [1]: (interreduced) ideal obtained by substituting from i those vars
       appearing in p, which occur in the linear part of i (or which can
       be expressed directly in the remaining variables, if e!=0)
  [2]: ideal, variables which have been substituted
  [3]: ideal, i-th element defines substitution of i-th var in [2]
  [4]: ideal of variables of basering, substituted ones are set to 0
  [5]: ideal, describing the map from the basering, say k[x(1..m)], to
       itself onto k[..variables fom [4]..] and [1] is the image of i
  [6]: int, # of vars considered for substitution (= # of factors of p)
  @end format
  The ideal i is generated by [1] and [3] in k[x(1..m)], the map [5]
  maps [3] to 0, hence induces an isomorphism
  @format
            k[x(1..m)]/i -> k[..variables fom [4]..]/[1]
  @end format
NOTE:    the proc uses @code{execute} to create a ring with ordering dp and vars
         placed correctly and then applies @code{elimpart}.
EXAMPLE: example elimpartanyr; shows an example
"
{
   def P = basering;
   int j,n,e = 0,0,1;
   poly p = product(maxideal(1));
   if ( size(#)==1 ) { p=#[1]; }
   if ( size(#)==2 ) { p=#[1]; e=#[2]; }
   string a,b;
   for ( j=1; j<=nvars(P); j++ )
   {
      if (deg(p/var(j))>=0) { a = a+varstr(j)+","; n = n+1; }
      else { b = b+varstr(j)+","; }
   }
   if ( size(b) != 0 ) { b = b[1,size(b)-1]; }
   else { a = a[1,size(a)-1]; }
   execute("ring gnir ="+charstr(P)+",("+a+b+"),dp;");
   ideal i = imap(P,i);
   list L = elimpart(i,n,e)+list(n);
   setring P;
   list L = imap(gnir,L);
   return(L);
}
example
{ "EXAMPLE:"; echo = 2;
   ring s=0,(x,y,z),dp;
   ideal i = x3+y2+z,x2y2+z3,y+z+1;
   elimpartanyr(i,z);
}
///////////////////////////////////////////////////////////////////////////////

proc fastelim (ideal i, poly p, list #)
"USAGE:   fastelim(i,p[h,o,a,b,e,m]); i=ideal, p=polynomial; h,o,a,b,e=integers
          p: product of variables to be eliminated;@*
  Optional parameters:
  @format
  - h !=0: use Hilbert-series driven std-basis computation
  - o !=0: use proc @code{valvars} for a - hopefully - optimal ordering of vars
  - a !=0: order vars to be eliminated w.r.t. increasing complexity
  - b !=0: order vars not to be eliminated w.r.t. increasing complexity
  - e !=0: use @code{elimpart} first to eliminate easy part
  - m !=0: compute a minimal system of generators
  @end format
  (default: h,o,a,b,e,m = 0,1,0,0,0,0)
RETURN:  ideal obtained from i by eliminating those variables, which occur in p
EXAMPLE: example fastelim; shows an example.
"
{
   def P = basering;
   int h,o,a,b,e,m = 0,1,0,0,0,0;
   if ( size(#) == 1 ) { h=#[1]; }
   if ( size(#) == 2 ) { h=#[1]; o=#[2]; }
   if ( size(#) == 3 ) { h=#[1]; o=#[2]; a=#[3]; }
   if ( size(#) == 4 ) { h=#[1]; o=#[2]; a=#[3]; b=#[4];}
   if ( size(#) == 5 ) { h=#[1]; o=#[2]; a=#[3]; b=#[4]; e=#[5]; }
   if ( size(#) == 6 ) { h=#[1]; o=#[2]; a=#[3]; b=#[4]; e=#[5]; m=#[6]; }
   list L = elimpartanyr(i,p,e);
   poly q = product(L[2]);     //product of vars which are already eliminated
   if ( q==0 ) { q=1; }
   p = p/q;                    //product of vars which must still be eliminated
   int nu = size(L[5])-size(L[2]);   //number of vars which must still be eliminated
   if ( p==1 )                 //ready if no vars are left
   {                           //compute minbase if 3-rd argument !=0
      if ( m != 0 ) { L[1]=minbase(L[1]); }
      return(L);
   }
//---------------- create new ring with remaining variables -------------------
   string newvar = string(L[4]);
   L = L[1],p;
   execute("ring r1=("+charstr(P)+"),("+newvar+"),"+"dp;");
   list L = imap(P,L);
//------------------- find "best" ordering of variables  ----------------------
   newvar = string(maxideal(1));
   if ( o != 0 )
   {
      list ordevar = valvars(L[1],a,L[2],b);
      intvec v = ordevar[1];
      newvar=string(sort(maxideal(1),v)[1]);
//------------ create new ring with "best" ordering of variables --------------
      changevar("r0",newvar);
      list L = imap(r1,L);
      kill r1;
      def r1 = r0;
      kill r0;
   }
//----------------- h==0: eliminate remaining vars directly -------------------
   if ( h == 0 )
   {
      L[1] = eliminate(L[1],L[2]);
      def r2 = r1;
   }
   else
//------- h!=0: homogenize and compute Hilbert-series using hilbvec ----------
   {
      intvec hi = hilbvec(L[1]);         // Hilbert-series of i
      execute("ring r2=("+charstr(P)+"),("+varstr(basering)+",@homo),dp;");
      list L = imap(r1,L);
      L[1] = homog(L[1],@homo);          // @homo = homogenizing var
//---- use Hilbert-series to eliminate variables with Hilbert-driven std -----
      L[1] = eliminate(L[1],L[2],hi);
      L[1]=subst(L[1],@homo,1);          // dehomogenize by setting @homo=1
   }
   if ( m != 0 )                         // compute minbase
   {
      if ( #[1] != 0 ) { L[1] = minbase(L[1]); }
   }
   def id = L[1];
   setring P;
   return(imap(r2,id));
}
example
{ "EXAMPLE:"; echo = 2;
   ring s=31991,(e,f,x,y,z,t,u,v,w,a,b,c,d),dp;
   ideal i = w2+f2-1, x2+t2+a2-1,  y2+u2+b2-1, z2+v2+c2-1,
            d2+e2-1, f4+2u, wa+tf, xy+tu+ab;
   fastelim(i,xytua,1,1);       //with hilb,valvars
   fastelim(i,xytua,1,0,1);     //with hilb,minbase
}
///////////////////////////////////////////////////////////////////////////////

proc faststd (@id,string @s1,string @s2, list #)
"USAGE:   faststd(id,s1,s2[,\"hilb\",\"sort\",\"dec\",o,\"blocks\"]);
         id=ideal/module, s1,s2=strings (names for new ring and maped id)
         o = string (allowed ordstring:\"lp\",\"dp\",\"Dp\",\"ls\",\"ds\",\"Ds\")
         \"hilb\",\"sort\",\"dec\",\"block\" options for Hilbert-std, sortandmap
COMPUTE: create a new ring (with \"best\" ordering of vars) and compute a
         std-basis of id (hopefully faster)
         - If say, s1=\"R\" and s2=\"j\", the new basering has name R and the
           std-basis of the image of id in R has name j
         - \"hilb\"  : use Hilbert-series driven std-basis computation
         - \"sort\"  : use 'sortandmap' for a best ordering of vars
         - \"dec\"   : order vars w.r.t. decreasing complexity (with \"sort\")
         - \"block\" : create blockordering, each block having ordstr=o, s.t.
                     vars of same complexity are in one block (with \"sort\")
         - o       : defines the basic ordering of the resulting ring
         default: o = ordering of 1-st block of basering - if it is allowed,
                      else o=\"dp\",
                  \"sort\", if none of the optional parameters is given
RETURN:  nothing
NOTE:    This proc is only useful for hard problems where other methods fail.
         \"hilb\" is useful for hard orderings (as \"lp\") or for characteristic 0,
         it is correct for \"lp\",\"dp\",\"Dp\" (and for blockorderings combining
         these) but not for s-orderings or if the vars have different weights.
         There seem to be only few cases in which \"dec\" is fast
EXAMPLE: example faststd; shows an example.
"
{
   def @P = basering;
   int @h,@s,@n,@m,@ii = 0,0,0,0,0;
   string @o,@va,@c = ordstr(basering),"","";
//-------------------- prepare ordering and set options -----------------------
   if ( @o[1]=="c" or @o[1]=="C")
      {  @o = @o[3,2]; }
   else
      { @o = @o[1,2]; }
   if( @o[1]!="d" and @o[1]!="D" and @o[1]!="l")
      { @o="dp"; }

   if (size(#) == 0 )
      { @s = 1; }
   for ( @ii=1; @ii<=size(#); @ii++ )
   {
      if ( typeof(#[@ii]) != "string" )
      {
         "// wrong syntax! type: help faststd";
         return();
      }
      else
      {
         if ( #[@ii] == "hilb"  ) { @h = 1; }
         if ( #[@ii] == "dec"   ) { @n = 1; }
         if ( #[@ii] == "block" ) { @m = 1; }
         if ( #[@ii] == "sort"  ) { @s = 1; }
         if ( #[@ii]=="lp" or #[@ii]=="dp" or #[@ii]=="Dp" or #[@ii]=="ls"
              or #[@ii]=="ds" or #[@ii]=="Ds" ) { @o = #[@ii]; }
      }
   }
   if( voice==2 ) { "// choosen options, hilb sort dec block:",@h,@s,@n,@m; }

//-------------------- nosort: create ring with new name ----------------------
   if ( @s==0 )
   {
      execute("ring "+@s1+"=("+charstr(@P)+"),("+varstr(@P)+"),("+@o+");");
      def @id = imap(@P,@id);
      verbose(noredefine);
      def @P = basering;
      verbose(redefine);
      kill `@s1`;
      if ( @h==0 ) { @id = std(@id); }
   }
//--------- sort: create new ring with "best" ordering of variables -----------
 proc bestorder(@id,string @s1,string @s2,int @n,string @o,int @m,int @l)
  {
     intvec getoption = option(get);
     option(redSB);
     @id = interred(sort(@id)[1]);
     poly @p = product(maxideal(1),1..@l);
     def i,s1,s2,n,p,o,m = @id,@s1,@s2,@n,@p,@o,@m;
     sortandmap(i,s1,s2,n,p,0,o,m);
     option(set,getoption);
     keepring(basering);
  }
//---------------------- no hilb: compute SB directly -------------------------
   if ( @s != 0 and @h == 0 )
   {
      bestorder(@id,@s1,@s2,@n,@o,@m,nvars(@P));
      verbose(noredefine);
      def @P = basering;
      verbose(redefine);
      kill `@s1`;
      def @id = `@s2`;
      @id = std(@id);
   }
//------- hilb: homogenize and compute Hilbert-series using hilbvec -----------
// this uses another standardbasis computation
   if ( @h != 0 )
   {
      execute("ring @Q=("+charstr(@P)+"),("+varstr(@P)+",@homo),("+@o+");");
      def @id = imap(@P,@id);
      kill @P;
      @id = homog(@id,@homo);               // @homo = homogenizing var
      if ( @s != 0 )
      {
         bestorder(@id,@s1,@s2,@n,@o,@m,nvars(@Q)-1);
         verbose(noredefine);
         def @Q= basering;
         kill `@s1`;
         def @id = `@s2`;
         verbose(redefine);
      }
      intvec @hi;                     // encoding of Hilbert-series of i
      @hi = hilbvec(@id);
      //if ( @s!=0 ) { @hi = hilbvec(@id,"32003",ordstr(@Q)); }
      //else { @hi = hilbvec(@id); }
//-------------------------- use Hilbert-driven std --------------------------
      @id = std(@id,@hi);
      @id = subst(@id,@homo,1);             // dehomogenize by setting @homo=1
      @va = varstr(@Q)[1,size(varstr(@Q))-6];
      if ( @s!=0 )
      {
         @o = ordstr(@Q);
         if ( @o[1]=="c" or @o[1]=="C") { @o = @o[1,size(@o)-6]; }
         else { @o = @o[1,size(@o)-8] + @o[size(@o)-1,2]; }
      }
      execute("ring @P=("+charstr(@Q)+"),("+@va+"),("+@o+");");
      def @id = imap(@Q,@id);
   }
   def `@s1` = @P;
   def `@s2` = @id;
   attrib(`@s2`,"isSB",1);
   export(`@s2`);
   kill @P;
   keepring(basering);
   if( voice==2 ) { "// basering is now "+@s1+", std-basis has name "+@s2; }
   return();
}
example
{ "EXAMPLE:"; echo = 2;
   ring s = 0,(e,f,x,y,z,t,u,v,w,a,b,c,d),(c,lp);
   ideal i = w2+f2-1, x2+t2+a2-1,  y2+u2+b2-1, z2+v2+c2-1,
            d2+e2-1, f4+2u, wa+tf, xy+tu+ab;
  option(prot); timer=1;
   int time = timer;
   ideal j=std(i);
   timer-time;
   show(R);dim(j),mult(j);
   int time = timer;
   faststd(i,"R","i");                      // use "best" ordering of vars
   timer-time;
   show(R);dim(i),mult(i);
   setring s;time = timer;
   faststd(i,"R","i","hilb");                // hilb-std only
   timer-time;
   show(R);dim(i),mult(i);
   setring s;time = timer;
   faststd(i,"R","i","hilb","sort");         // hilb-std,"best" ordering
   timer-time;
   show(R);dim(i),mult(i);
    setring s;time = timer;
   faststd(i,"R","i","hilb","sort","block","dec"); // hilb-std,"best",blocks
   timer-time;
   show(R);dim(i),mult(i);
  setring s;time = timer;
   timer-time;time = timer;
   faststd(i,"R","i","sort","block","Dp"); //"best",decreasing,Dp-blocks
   timer-time;
   show(R);dim(i),mult(i);
}
///////////////////////////////////////////////////////////////////////////////

proc findvars(id, list #)
"USAGE:   findvars(id [,any] ); id=poly/ideal/vector/module/matrix, any=any type
RETURN:  if no second argument is present: ideal of variables occuring in id,@*
         if a second argument is given (of any type): list L with 4 entries:
  @format
  L[1]: ideal of variables occuring in id
  L[2]: intvec of variables occuring in id
  L[3]: ideal of variables not occuring in id
  L[4]: intvec of variables not occuring in id
  @end format
EXAMPLE: example findvars; shows an example
"
{
   int ii,n;
   ideal found, notfound;
   intvec f,nf;
   n = nvars(basering);
   ideal i = simplify(ideal(matrix(id)),10);
   matrix M[ncols(i)][1] = i;
   vector v = module(M)[1];
   ideal max = maxideal(1);

   for (ii=1; ii<=n; ii++)
   {
      if ( v != subst(v,var(ii),0) )
      {
         found = found+var(ii);
         f = f,ii;
      }
      else
      {
         notfound = notfound+var(ii);
         nf = nf,ii;
      }
   }
   if ( size(f)>1 ) { f = f[2..size(f)]; }      //intvec of found vars
   if ( size(nf)>1 ) { nf = nf[2..size(nf)]; }  //intvec of vars not found
   if( size(#)==0 )  { return(found); }
   if( size(#)!=0 )  { list L = found,f,notfound,nf; return(L); }
}
example
{ "EXAMPLE:"; echo = 2;
   ring s  = 0,(e,f,x,y,t,u,v,w,a,d),dp;
   ideal i = w2+f2-1, x2+t2+a2-1;
   findvars(i);
   findvars(i,1);
}
///////////////////////////////////////////////////////////////////////////////

proc hilbvec (@id, list #)
"USAGE:   hilbvec(id[,c,o]); id=poly/ideal/vector/module/matrix, c,o=strings,@*
          c=char, o=ordering used by @code{hilb}@*
          (default: c=\"32003\", o=\"dp\")
RETURN:  intvec of 1-st Hilbert-series of id, computed in char c and ordering o
NOTE:    id must be homogeneous (i.e. all vars have weight 1)
EXAMPLE: example hilbvec; shows an example
"
{
   def @P = basering;
   string @c,@o = "32003", "dp";
   if ( size(#) == 1 ) {  @c = #[1]; }
   if ( size(#) == 2 ) {  @c = #[1]; @o = #[2]; }
   string @si = typeof(@id)+" @i = "+string(@id)+";";  //** weg
   execute("ring @r=("+@c+"),("+varstr(basering)+"),("+@o+");");
   //**def i = imap(P,@id);
   execute(@si);                   //** weg
   //show(basering);
   @i = std(@i);
   intvec @hi = hilb(@i,1);         // intvec of 1-st Hilbert-series of id
   return(@hi);
}
example
{ "EXAMPLE:"; echo = 2;
   ring s   = 0,(e,f,x,y,z,t,u,v,w,a,b,c,d,H),dp;
   ideal id = w2+f2-1, x2+t2+a2-1,  y2+u2+b2-1, z2+v2+c2-1,
              d2+e2-1, f4+2u, wa+tf, xy+tu+ab;
   id = homog(id,H);
   hilbvec(id);
}
///////////////////////////////////////////////////////////////////////////////

proc linearpart (id)
"USAGE:   linearpart(id);  id=ideal/module
RETURN:  generators of id of total degree <= 1
EXAMPLE: example linearpart; shows an example
"
{
   return(degreepart(id,0,1));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   ideal i=1+x+x2+x3,3,x+3y+5z;
   linearpart(i);
   module m=[x,y,z],x*[x3,y2,z],[1,x2,z3,0,1];
   show(linearpart(m));
}
///////////////////////////////////////////////////////////////////////////////

proc tolessvars (id ,list #)
"USAGE:   tolessvars(id [,s1,s2] ); id poly/ideal/vector/module/matrix,
          s1,s2=strings@*
          s1: name of new ring,@*
          s2: new ordering@*
          (default: s1=\"R(n)\" where n is the # of vars in the new ring,
          s2=\"dp\" or \"ds\" depending whether the first block of the old
          ordering is a p- resp. an s-ordering)

CREATE:  nothing, if id contains all vars of the basering.@*
         Else, create a ring with same char as the basering, but possibly less
         variables (only those variables which actually occur in id) and map
         id to the new ring, which will be the basering after the proc has
         finished.
DISPLAY: If printlevel >=0, display ideal of vars, which have been ommitted from
         the old ring
RETURN:  the original ideal id (see NOTE)
NOTE:    You must not type, say, 'ideal id=tolessvars(id);' since the ring
         to which 'id' would belong will only be defined by the r.h.s.. But you
         may type 'def id=tolessvars(id);' or 'list id=tolessvars(id);'
         since then 'id' does not a priory belong to a ring, its type will
         be defined by the right hand side. Moreover, do not use a name which
         occurs in the old ring, for the same reason.
EXAMPLE: example tolessvars; shows an example
"
{
//---------------- initialisation and check occurence of vars -----------------
   int s,ii,n,fp,fs;
   string s1,s2,newvar;
   int pr = printlevel-voice+3;  // p = printlevel+1 (default: p=1)
   def P = basering;
   s2 = ordstr(P);

   list L = findvars(id,1);
   newvar = string(L[1]);    // string of new variables
   n = size(L[1]);           // number of new variables
   if( n == 0 )
   {
      dbprint( pr,"","// no variable occured in "+typeof(id)+", no change of ring!");
      return(id);
   }
   if( n == nvars(P) )
   {
      dbprint( pr,"","// all variables occured in "+typeof(id)+", no change of ring!");
      return(id);
   }
//----------------- prepare new ring, map to it and return --------------------
   s1 = "R("+string(n)+")";
   if ( size(#) == 0 )
   {
       fp = find(s2,"p");
       fs = find(s2,"s");
       if( fs==0 or (fs>=fp && fp!=0) ) { s2="dp"; }
       else {  s2="ds"; }
   }
   if ( size(#) ==1 ) { s1=#[1]; }
   if ( size(#) ==2 ) { s1=#[1]; s2=#[2]; }
   //dbprint( pr,"","// variables which did not occur:",simplify(max,2) );
   dbprint( pr,"","// variables which did not occur:",L[3] );

   execute("ring "+s1+"=("+charstr(P)+"),("+newvar+"),("+s2+");");
   def id = imap(P,id);
   export(basering);
   keepring (basering);
   dbprint( pr,"// basering is now "+s1 );
   return(id);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r  = 0,(x,y,z),dp;
   ideal i = y2-x3,x-3,y-2x;
   def j   = tolessvars(i,"R_r","lp");
   show(basering);
   j;
   kill R_r;
}
///////////////////////////////////////////////////////////////////////////////

proc solvelinearpart (id,list #)
"USAGE:   solvelinearpart(id [,n] );  id=ideal/module, n=integer,@*
          (default: n=0)
RETURN:  (interreduced) generators of id of degree <=1 in reduced triangular
         form if n=0 [non-reduced triangular form if n!=0]
ASSUME:  monomial ordering is a global ordering (p-ordering)
NOTE:    may be used to solve a system of linear equations
         see proc @code{gauss_row} from 'matrix.lib' for a different method
WARNING: the result is very likely to be false for 'real' coefficients, use
         char 0 instead!
EXAMPLE: example solvelinearpart; shows an example
"
{
   intvec getoption = option(get);
   option(redSB);
   if ( size(#)!=0 )
   {
      if(#[1]!=0) { option(noredSB); }
   }
   def lin = interred(degreepart(id,0,1));
   if ( size(#)!=0 )
   {
      if(#[1]!=0)
      {
         return(lin);
      }
   }
   option(set,getoption);
   return(simplify(lin,1));
}
example
{ "EXAMPLE:"; echo = 2;
   // Solve the system of linear equations:
   //         3x +   y +  z -  u = 2
   //         3x +  8y + 6z - 7u = 1
   //        14x + 10y + 6z - 7u = 0
   //         7x +  4y + 3z - 3u = 3
   ring r = 0,(x,y,z,u),lp;
   ideal i= 3x +   y +  z -  u,
           13x +  8y + 6z - 7u,
           14x + 10y + 6z - 7u,
            7x +  4y + 3z - 3u;
   ideal j= 2,1,0,3;
   j = i-j;                        // difference of 1x4 matrices
                                   // compute reduced triangular form, setting
   solvelinearpart(j);             // the RHS equal 0 gives the solutions!
   solvelinearpart(j,1); "";       // triangular form, not reduced
}
///////////////////////////////////////////////////////////////////////////////

proc sortandmap (@id,string @s1,string @s2, list #)
"USAGE:   sortandmap(id,s1,s2[,n1,p1,n2,p2...,o1,m1,o2,m2...]);@*
         id=poly/ideal/vector/module,@*
         s1,s2 = strings (names for new ring and mapped id),@*
         p1,p2,...= polynomials (product of variables),@*
         n1,n2,...= integers,@*
         o1,o2,...= strings,@*
         m1,m2,...= integers@*
         (default: p1=product of all vars, n1=0, o1=\"dp\",m1=0)
         the last pi (containing the remaining vars) may be omitted
CREATE:  a new ring and map id into it, the new ring has same char as basering
         but with new ordering and vars sorted in the following manner:
  @format
  - each block of vars occuring in pi is sorted w.r.t. its complexity in id,
  - ni controls the sorting in i-th block (= vars occuring in pi):
    ni=0 (resp.!=0) means that less (resp. more) complex vars come first
  - oi and mi define the monomial ordering of the i-th block:
    if mi =0, oi=ordstr(i-th block)
    if mi!=0, the ordering of the i-th block itself is a blockordering,
      each subblock having ordstr=oi, such that vars of same complexity are
      in one block
  @end format
  Note that only simple ordstrings oi are allowed:
  \"lp\",\"dp\",\"Dp\",\"ls\",\"ds\",\"Ds\".
RETURN:  nothing
NOTE:    We define a variable x to be more complex than y (with respect to id)
         if val(x) > val(y) lexicographically, where val(x) denotes the
         valuation vector of x:@*
         consider id as list of polynomials in x with coefficients in the
         remaining variables. Then:@*
         val(x) = (maximal occuring power of x,  # of all monomials in leading
         coefficient, # of all monomials in coefficient of next smaller power
         of x,...).
EXAMPLE: example sortandmap; shows an example
"
{
   def @P = basering;
   int @ii,@jj;
   intvec @v;
   string @o;
//----------------- find o in # and split # into 2 lists ---------------------
   # = # +list("dp",0);
   for ( @ii=1; @ii<=size(#); @ii++)
   {
      if ( typeof(#[@ii])=="string" )  break;
   }
   if ( @ii==1 ) { list @L1 = list(); }
   else { list @L1 = #[1..@ii-1]; }
   list @L2 = #[@ii..size(#)];
   list @L = sortvars(@id,@L1);
   string @va = string(@L[1]);
   list @l = @L[2];   //e.g. @l[4]=intvec describing permutation of 1-st block
//----------------- construct correct ordering with oi and mi ----------------
   for ( @ii=4; @ii<=size(@l); @ii=@ii+4 )
   {
      @L2=@L2+list("dp",0);
      if ( @L2[@ii/2] != 0)
      {
         @v = @l[@ii];
         for ( @jj=1; @jj<=size(@v); @jj++ )
         {
           @o = @o+@L2[@ii/2 -1]+"("+string(@v[@jj])+"),";
         }
      }
      else
      {
         @o = @o+@L2[@ii/2 -1]+"("+string(size(@l[@ii/2]))+"),";
      }
   }
   @o=@o[1..size(@o)-1];
//------------------ create new ring and make objects global -----------------
   execute("ring "+@s1+"=("+charstr(@P)+"),("+@va+"),("+@o+");");
   def @id = imap(@P,@id);
   execute("def "+ @s2+"=@id;");
   execute("export("+@s1+");");
   execute("export("+@s2+");");
   keepring(basering);
   return();
}
example
{ "EXAMPLE:"; echo = 2;
   ring s = 32003,(x,y,z),dp;
   ideal i=x3+y2,xz+z2;
   sortandmap(i,"R_r","i");
   // i is now an ideal in the new basering R_r
   show(R_r);
   kill R_r; setring s;
   sortandmap(i,"R_r","i",1,xy,0,z,0,"ds",0,"lp",0);
   show(R_r);
   kill R_r;
}
///////////////////////////////////////////////////////////////////////////////

proc sortvars (id, list #)
"USAGE:   sortvars(id[,n1,p1,n2,p2,...]);@*
         id=poly/ideal/vector/module,@*
         p1,p2,...= polynomials (product of vars),@*
         n1,n2,...=integers@*
         (default: p1=product of all vars, n1=0)
         the last pi (containing the remaining vars) may be omitted
COMPUTE: sort variables with respect to their complexity in id
RETURN:  list of two elements, an ideal and a list:
  @format
  [1]: ideal, variables of basering sorted w.r.t their complexity in id
       ni controls the ordering in i-th block (= vars occuring in pi):
       ni=0 (resp.!=0) means that less (resp. more) complex vars come first
  [2]: a list with 4 entries for each pi:
       ideal ai : vars of pi in correct order,
       intvec vi: permutation vector describing the ordering in ai,
       intmat Mi: valuation matrix of ai, the columns of Mi being the
                  valuation vectors of the vars in ai
       intvec wi: size of 1-st, 2-nd,... block of identical columns of Mi
                  (vars with same valuation)
  @end format
NOTE:    We define a variable x to be more complex than y (with respect to id)
         if val(x) > val(y) lexicographically, where val(x) denotes the
         valuation vector of x:@*
         consider id as list of polynomials in x with coefficients in the
         remaining variables. Then:@*
         val(x) = (maximal occuring power of x,  # of all monomials in leading
         coefficient, # of all monomials in coefficient of next smaller power
         of x,...).
EXAMPLE: example sortvars; shows an example
"
{
   int ii,jj,n,s;
   list L = valvars(id,#);
   list L2, L3 = L[2], L[3];
   list K; intmat M; intvec v1,v2,w;
   ideal i = sort(maxideal(1),L[1])[1];
   for ( ii=1; ii<=size(L2); ii++ )
   {
      M = transpose(L3[2*ii]);
      M = M[L2[ii],1..nrows(L3[2*ii])];
      w = 0; s = 0;
      for ( jj=1; jj<=nrows(M)-1; jj++ )
      {
         v1 = M[jj,1..ncols(M)];
         v2 = M[jj+1,1..ncols(M)];
         if ( v1 != v2 ) { n=jj-s; s=s+n; w = w,n; }
      }
      w=w,nrows(M)-s; w=w[2..size(w)];
      K = K+sort(L3[2*ii-1],L2[ii])+list(transpose(M))+list(w);
   }
   L = i,K;
   return(L);
}
example
{ "EXAMPLE:"; echo = 2;
   ring s=0,(x,y,z,w),dp;
   ideal i = x3+y2+yw2,xz+z2,xyz-w2;
   sortvars(i,0,xy,1,zw);
}
///////////////////////////////////////////////////////////////////////////////

proc valvars (id, list #)
"USAGE:   valvars(id[,n1,p1,n2,p2,...]);@*
         id=poly/ideal/vector/module,@*
         p1,p2,...= polynomials (product of vars),@*
         n1,n2,...= integers,

         ni controls the ordering of vars occuring in pi: ni=0 (resp.!=0) means
         that less (resp. more) complex vars come first@*
         (default: p1=product of all vars, n1=0)
         the last pi (containing the remaining vars) may be omitted
COMPUTE: valuation (complexity) of variables with respect to id.@*
         ni controls the ordering of vars occuring in pi:@*
         ni=0 (resp.!=0) means that less (resp. more) complex vars come first.
RETURN:  list with 3 entries:
  @format
  [1]: intvec, say v, describing the permutation such that the permuted
       ringvariables are ordered with respect to their complexity in id
  [2]: list of intvecs, i-th intvec, say v(i) describing permutation
       of vars in a(i) such that v=v(1),v(2),...
  [3]: list of ideals and intmat's, say a(i) and M(i), where
       a(i): factors of pi,
       M(i): valuation matrix of a(i), such that the j-th column of M(i)
             is the valuation vector of j-th generator of a(i)
         @end format
NOTE:    Use @code{sortvars} in order to actually sort the variables!
         We define a variable x to be more complex than y (with respect to id)
         if val(x) > val(y) lexicographically, where val(x) denotes the
         valuation vector of x:@*
         consider id as list of polynomials in x with coefficients in the
         remaining variables. Then:@*
         val(x) = (maximal occuring power of x,  # of all monomials in leading
         coefficient, # of all monomials in coefficient of next smaller power
         of x,...).
EXAMPLE: example valvars; shows an example
"
{
//---------------------------- initialization ---------------------------------
   int ii,jj,kk,n;
   list L;                    // list of valuation vectors in one block
   intvec vec;                // describes permutation of vars (in one block)
   list blockvec;             // i-th element = vec of i-th block
   intvec varvec;             // result intvector
   list Li;                   // result list of ideals
   list LM;                   // result list of intmat's
   intvec v,w,s;              // w valuation vector for one variable
   matrix C;                  // coefficient matrix for different variables
   ideal i = simplify(ideal(matrix(id)),10);

//---- for each pii in # create ideal a(ii) intvec v(ii) and list L(ii) -------
// a(ii) = ideal of vars in product, v(ii)[j]=k <=> a(ii)[j]=var(k)

   v = 1..nvars(basering);
   int l = size(#);
   if ( l >= 2 )
   {
      ideal m=maxideal(1);
      for ( ii=2; ii<=l; ii=ii+2 )
      {
         int n(ii) = #[ii-1];
         ideal a(ii);
         intvec v(ii);
         for ( jj=1; jj<=nvars(basering); jj++ )
         {
            if ( #[ii]/var(jj) != 0)
            {
               a(ii) = a(ii) + var(jj);
               v(ii)=v(ii),jj;
               m[jj]=0;
               v[jj]=0;
            }
         }
         v(ii)=v(ii)[2..size(v(ii))];
      }
      if ( size(m)!=0 )
      {
         l = 2*(l/2)+2;
         ideal a(l) = simplify(m,2);
         intvec v(l) = compress(v);
         int n(l);
         if ( size(#)==l-1 ) { n(l) = #[l-1]; }
      }
   }
   else
   {
      l = 2;
      ideal a(2) = maxideal(1);
      intvec v(2) = v;
      int n(2);
      if ( size(#)==1 ) { n(2) = #[1]; }
   }
//------------- start loop to order variables in each a(ii) -------------------

   for ( kk=2; kk<=l; kk=kk+2 )
   {
      L = list();
      n = 0;
//---------------- get valuation of all variables in a(kk) --------------------
      for ( ii=1; ii<=size(a(kk)); ii++ )
      {
         C = coeffs(i,a(kk)[ii]);
         w = nrows(C); // =(maximal occuring power of a(kk)[ii])+1
         for ( jj=w[1]; jj>1; jj-- )
         {
            s = size(C[jj,1..ncols(C)]);
            w[w[1]-jj+2] = sum(s);
         }
         // w[1] should represent the maximal occuring power of a(kk)[ii] so it
         // has to be decreased by 1 since otherwise the constant term is also
         // counted
         w[1]=w[1]-1;

         L[ii]=w;
         n = size(w)*(size(w) > n) + n*(size(w) <= n);
      }
      intmat M(kk)[size(a(kk))][n];
      for ( ii=1; ii<=size(a(kk)); ii++ )
      {
         if ( n==1 ) { w = L[ii]; M(kk)[ii,1] = w[1]; }
         else  { M(kk)[ii,1..n] = L[ii]; }
      }
      LM[kk-1] = a(kk);
      LM[kk] = transpose(compress(M(kk)));
//------------------- compare valuation and insert in vec ---------------------
      vec = sort(L)[2];
      if ( n(kk) != 0 ) { vec = vec[size(vec)..1]; }
      blockvec[kk/2] = vec;
      vec = sort(v(kk),vec)[1];
      varvec = varvec,vec;
   }
   varvec = varvec[2..size(varvec)];
   list result = varvec,blockvec,LM;
   return(result);
}
example
{ "EXAMPLE:"; echo = 2;
   ring s=0,(x,y,z,a,b),dp;
   ideal i=ax2+ay3-b2x,abz+by2;
   valvars (i,0,xyz);
}
///////////////////////////////////////////////////////////////////////////////
/*

 ring s=31991,(e,f,x,y,z,t,u,v,w,a,b,c,d),dp;
 ring s=31991,(x,y,z,t,u,v,w,a,b,c,d,f,e,h),dp; //standard
 ring s1=31991,(y,u,b,c,a,z,t,x,v,d,w,e,f,h),dp; //gut
v;
13,12,11,10,8,7,6,5,4,3,2,1,9,14
print(matrix(sort(maxideal(1),v)));
f,e,w,d,x,t,z,a,c,b,u,y,v,h
print(matrix(maxideal(1)));
y,u,b,c,a,z,t,x,v,d,w,e,f,h
v0;
14,9,12,11,10,8,7,6,5,4,3,2,1,13
print(matrix(sort(maxideal(1),v0)));
h,v,e,w,d,x,t,z,a,c,b,u,y,f
v1;v2;
9,12,11,10,8,7,6,5,4,3,2,1,13,14
13,12,11,10,8,7,6,5,4,3,2,1,9,14

Ev. Gute Ordnung fuer i:
========================
i=ad*x^d+ad-1*x^(d-1)+...+a1*x+a0, ad!=0
mit ar=(ar1,...,ark), k=size(i)
    arj in K[..x^..]
d=deg_x(i) := max{deg_x(i[k]) | k=1..size(i)}
size_x(i,deg_x(i)..0) := size(ad),...,size(a0)
x>y  <==
  1. deg_x(i)>deg_y(i)
  2. "=" in 1. und size_x lexikographisch

hier im Beispiel:
f: 5,1,0,1,2

u: 3,1,4

y: 3,1,3
b: 3,1,3
c: 3,1,3
a: 3,1,3
z: 3,1,3
t: 3,1,3

x: 3,1,2
v: 3,1,2
d: 3,1,2
w: 3,1,2
e: 3,1,2
probier mal:
 ring s=31991,(f,u,y,z,t,a,b,c,v,w,d,e,h),dp; //standard

*/