//
version="$Id: fastsolv.lib,v 1.5 1998-05-29 15:01:56 Singular Exp $";
info="";

// ====================================================================
//  library fast.lib
//
// ====================================================================

proc Tsimplify
"USAGE:    simplify(id,n);  id ideal, n integer
RETURNS:  list of two ideals the first containing the simplified
          elements, the second the variables which already have
          been eliminated
ASSUME:   basering has ordering dp, the elimination of the variables
          n+1,... is not supported
EXAMPLE:  example simplify; shows an example
"
{
  def bsr=basering;
  int @n, @m, @k;
  option (redSB);
  ideal @i = interred(#[1]);
  ideal @j, @im, @va;
  for (@n=1;@n<=size(@i);@n=@n+1)
  {
    if (deg(@i[@n])==1)
    {
      @k=0;
      for (@m=#[2]+1;@m<=nvars(basering);@m=@m+1)
      {
        if (lead(@i[@n])/var(@m)!=0)
        {
          @k=1;
          break;
        }
      }
      if (@k==0)
      {
        @va=@va+lead(@i[@n]);
        @i[@n]=0;
      }
    }
  }
  @i=@i+@j;
  map @phi;
  for (@m=1;@m<=#[2];@m=@m+1)
  {
    for (@n=size(@i);@n>=1;@n=@n-1)
    {
      if (deg(@i[@n]/var(@m))==0)
      {
        @im = maxideal(1);
        @im[@m]=(-1/coeffs(@i[@n],var(@m))[2,1])*(@i[@n]-coeffs(@i[@n],var(@m))[2,1]*var(@m));
        @phi = basering,@im;
        @i=@phi(@i);
        @i=@i+@j;
        @va=@va+var(@m);
        break;// hier muesste evt. sorgfaeltiger ausgewaehlt werden
      }
    }
  }
  @i=interred(@i);
  option(noredSB);
  return(@i,@va);
}
example
{
   ring s=181,(x,y,z,t,u,v,w,a,b,c,d,f,e),dp;
   ideal j=
            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,
            yz+uv+bc,
            cd+ze,
            x+y+z+e+1,
            t+u+v+f-1,
            w+a+b+c+d;
    list l=Tsimplify(j,10);
    l;
}


proc simplifyInBadRings
"USAGE:    simplify(id,n);  id ideal, n product of the variables
          to be eliminated later
RETURNS:  list of two ideals the first containing the simplified
          elements, the second the variables which already have
          been eliminated
EXAMPLE:  example simplify; shows an example
"
{
  int @i,@j;
  string @es=",(";
  string @ns="),dp;";
  for (@i=1;@i<=nvars(basering);@i=@i+1)
  {
    if (deg(#[2]/var(@i))>=0)
    {
      @es=@es+varstr(@i)+",";
      @j=@j+1;
    }
    else
    {
      @ns=","+varstr(@i)+@ns;
    }
  }
  string @s=@es[1,size(@es)-1]+@ns;
  @ns=nameof(basering);
  ideal i=#[1];
  execute "ring @gnir ="+charstr(basering)+@s;
  ideal @laedi=imap(`@ns`,i);
  list @l=Tsimplify(@laedi,@j);
  ideal @l1, @l2=@l[1], @l[2];
  execute "setring "+@ns+";";
  return(imap(@gnir,@l1),imap(@gnir,@l2));
}
example
{
   ring s=0,(e,f,x,y,z,t,u,v,w,a,b,c,d),dp;
   ideal j=
            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,
            yz+uv+bc,
            cd+ze,
            x+y+z+e+1,
            t+u+v+f-1,
            w+a+b+c+d;
    list l=simplifyInBadRings(j,xyztuvwabc);
    l;
}


proc fastElim
"USAGE:    fastelim(id,v);  id ideal, v Product of variables to
          be eliminated.
RETURNS:  list of two ideals the first containing the simplified
          elements, the second the variables which already have
          been eliminated
EXAMPLE:  example fastElim; shows an example.
"{
   string @oldgnir = nameof(basering);  // Namen des alten Baserings merken
   string @varnames;
   int @i;
   poly @testpoly=1;      // Variablen, die schon eliminiert sind (nach Simplify)
   poly @elimvarspoly=#[2]; // Variablen, die noch eliminiert werden sollen

   // das Ideal vereinfachen. Beachte: l[1] erzeugt nicht mehr id!
   list @l=simplifyInBadRings(#[1],#[2]);
   ideal @a, @b = @l[1], @l[2];

   // ------------------------------------------------------------
   // in einen Ring abbilden, der die in (2) angegebenen Variablen
   // nicht mehr enthaelt.
   // ------------------------------------------------------------


   // Aus dem von simplifyInBadRings zurueckgegebenen Ideal(2), das die
   // Variablen anthealt, die eliminiert sind, ein Polynom bauen:
   for (@i=1; @i<=size(@b); @i=@i+1)
   {
     @testpoly=@testpoly*@b[@i];
     @elimvarspoly = subst(@elimvarspoly, @b[@i], 1);
   }

   //Wenn keine Variablen mehr uebrig sind, sind wir fertig:
   if (@elimvarspoly==1)
   {
     // Falls das 4. Argument 1 ist, eine minimale Basis berechnen
     if (#[4]==1)
     {
       @a=mres(@a,1)[1];
     }
     return(@a);
   }

   // Anhand dieses Polynoms ueberpruefen, welche Variablen unser
   // neuer Ring braucht: @varnames
   for (@i=1; @i<=nvars(basering); @i=@i+1)
   {
     if (@testpoly/var(@i)==0)
     {
       @varnames=@varnames + nameof(var(@i)) + ",";
     }
   }
   @varnames=@varnames[1, size(@varnames)-1];

   // Den neuen Ring erzeugen und das Ideal a darin abbilden:
   execute("ring @r1=" + charstr(basering) + ",(" + @varnames + ")," + "dp;" );
   ideal @a=imap(`@oldgnir`, @a);
   list #=imap(`@oldgnir`,#);

   // Nun die beste Reihenfolge fuer die Variablen berechnen:
   string @orderedvars=string(maxideal(1));
   if (#[3]==1)
   {
      poly @elimvarspoly=imap(`@oldgnir`,@elimvarspoly);
      intvec @v= valvars(@a,#[5], @elimvarspoly,#[6]);
      @v;
      @orderedvars=string(sort1(maxideal(1),@v));
   }

   // In einen Ring wechseln, der diese Reihenfolge beruecksichtig:
   execute("ring @r2=" + charstr(basering) + ",(" + @orderedvars + ",@homo),dp;" );
   ideal @a= imap(@r1, @a);
   kill @r1;

   //das ideal homogenisieren
   @a=homog(@a,@homo);

   // Die Variablen holen, die wirklich noch zu eliminieren sind, und
   // danach eliminate aufrufen.
   poly @elimvarspoly= imap(`@oldgnir`, @elimvarspoly);
   ideal @b=eliminate(@a, @elimvarspoly);

   // Falls das 4. Argument 1 ist, eine minimale Basis berechnen
   if (#[4]==1)
   {
     @b=mres(@b,1)[1];
   }

   // dehomogenisieren (@homo=1) und zurueck in den alten Ring
   @b=subst(@b,@homo,1);
   execute "setring "+@oldgnir+";";
   return(imap(@r2,@b));
}
example
{
   ring s=31991,(e,f,x,y,z,t,u,v,w,a,b,c,d),dp;
   ideal j=
            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,
            yz+uv+bc,
            cd+ze,
            x+y+z+e+1,
            t+u+v+f-1,
            w+a+b+c+d;
   fastElim(j,xyztuvwabcd,1,0,0,0);
}


proc sort1
{
  intvec @i=#[2];
  int @ii=size(@i);
  int @j;
  string @s=string(typeof(#[1]))+" @m;";
  execute(@s);
  for (@j=1;@j<=@ii;@j=@j+1)
  {
    @m[@j] = #[1][@i[@j]];
  }
  return(@m);
}


proc valvars
"
USAGE:    valvars(id, [n1, [m, [n2]]]); id (poly|ideal|vector|module);
                                        m monom; n1,n2 int
RETURNS:  an intvec describing the permutation such that the permuted
          ringvariables are ordered with respect to inreasing complexity
          [resp. decreasing complexity if n1 is present and non-zero].
          If m is present, the variables occuring in m are sorted in one
          block before all variables not occuring in m. n2 controls the
          order of this block as n1 does for the other one.
          A variable x is more complex than y (with respect to id) if, in the
          input, either the maximal occuring power of x is bigger than that
          of y or both are equal and the coefficient of x to this power
          contains more monomials than that of y.
EXAMPLE:  example valvars; shows an example
"
{
  int     @order1, @order2,     // order of vars not to elim resp. to elim
          @i1,                  // runs through all variables
          @i2,                  // # of vars to elim already found
          @N, @M,               // N = total # of vars, M = # of vars to elim
          @i, @j, @k, @size, @degree;
  intvec  @varvec, @degvec, @sizevec;
  poly    @m;
  ideal   @id;
  matrix  @C;

  @id = matrix( #[1] ); @m = 0;
  @N = nvars( basering ); @M = 0;
  @order1 = 1; @order2 = 1;
  if ( size(#) >= 4 )
    { @order2 = not( #[4] ) * 2 - 1; }
  if ( size(#) >= 3 )
    { @m = #[3]; @M = deg( #[3] ); }
  if ( size(#) >= 2 )
    { @order1 = not( #[2] ) * 2 - 1; }

  @i2 = 0;
  for ( @i1 = 1; @i1 <= @N; @i1++ )
  {
    // get valuation of current variable.
    @C = coeffs( @id, var( @i1 ));
    @degree = nrows( @C );
    @size = 0;
    for ( @j = ncols( @C ); @j >= 1; @j-- )
      { @size = @size + size( @C[ @degree, @j ] ); }

    // mind order and calculate limits for search
    if ( @M and size( @m / var( @i1 )))
    {
      @degree = @degree * @order2; @size = @size * @order2;
      @i2++;
      @i = @i2; @j = 1;
    }
    else
    {
      @degree = @degree * @order1; @size = @size * @order1;
      @i = @i1 - @i2 + @M; @j = @M+1;
    }

    // look where we have to insert the variable
    // "start = " + string( @j ) + ", end = " + string( @i );
    @degvec[ @i ] = @degree + 1; @sizevec[ @i ] = @size + 1;
    while ( @degvec[ @j ] < @degree )
      { @j++; }
    while ( @degvec[ @j ] == @degree and @sizevec[ @j ] < @size )
      { @j++; }

    // now shift the variables behind the current one back
    // "inserting at " + string( @j );
    for ( @k = @i-1; @k >= @j; @k-- )
    {
      @varvec[ @k+1 ] = @varvec[ @k ];
      @degvec[ @k+1 ] = @degvec[ @k ];
      @sizevec[ @k+1 ] = @sizevec[ @k ];
    }
    @varvec[ @j ] = @i1;             // print( @varvec );
    @degvec[ @j ] = @degree;         // print( @degvec );
    @sizevec[ @j ] = @size;          // print( @sizevec );
  }

  return( @varvec );
}
example
{
  echo = 1; ring @r=0,(a,b,c,x,y,z),lp;
  poly @f=a3+b4+xyz2+xyz+yz+1;
  valvars( @f );
  valvars( @f, 1 );
  valvars( @f, 0, abc );
  valvars( @f, 0, xyz );
  valvars( @f, 0, abc, 1 );
  valvars( @f, 1, abc, 1 );
  ideal @i=@f,z5,y5,-y5;
  valvars( @i );
  valvars( @i, 0, abcxyz );
  kill @r; echo = 0;
}