//////////////////////////////////////////////////////////////////////////////
version="$Id: sagbi.lib,v 1.13 2009-04-07 16:18:06 seelisch Exp $";
category="Commutative Algebra";
info="
LIBRARY:  sagbi.lib  Compute subalgebra bases analogous to Groebner bases for ideals
AUTHORS: Gerhard Pfister,        pfister@mathematik.uni-kl.de,
@*       Anen Lakhal,            alakhal@mathematik.uni-kl.de

PROCEDURES:
 sagbiRreduction(p,I);  perform one step subalgebra reducton (for short S-reduction) of p w.r.t I
 sagbiSPoly(I);   compute the S-polynomials of the Subalgebra defined by the genartors of I
 sagbiNF(id,I);   perform iterated S-reductions in order to compute Subalgebras normal forms
 sagbi(I);        construct SAGBI basis for the Subalgebra defined by I
 sagbiPart(I);    construct partial SAGBI basis for the Subalgebra defined by I
";

LIB "algebra.lib";
LIB "elim.lib";

///////////////////////////////////////////////////////////////////////////////
proc sagbiSPoly(id ,list #)
"USAGE: sagbiSPoly(id [,n]); id ideal, n positive integer.
RETURN: an ideal
@format
      - If (n=0 or default) an ideal, whose generators are the S-polynomials.
      - If (n=1) a list of size 2:
        the first element of this list is the ideal of S-polynomials.
        the second element of this list is the ring in which the ideal of algebraic
        relations is defined.
@end format
EXAMPLE: example sagbiSPoly; show an example "
{
  if(size(#)==0)
  {
    #[1]=0;
  }
  degBound=0;
  def bsr=basering;
  ideal vars=maxideal(1);
  ideal B=ideal(bsr);//when the basering is quotient ring this "type casting"
                    //gives th quotient ideal.
  int b=size(B);

  //In quotient rings,SINGULAR does not reduce polynomials w.r.t the
  //quotient ideal,therefore we should first 'reduce';if it is necessary for
  //computations to have a uniquely determined representant for each equivalent
  //class,which is the case of this procedure.

  if(b!=0)
  {
    id =reduce(id,groebner(0));
  }
  int n,m=nvars(bsr),ncols(id);
  int z;
  string mp=string(minpoly);
  ideal P;
  list L;

  if(id==0)
  {
    if(#[1]==0)
    {
      return(P);
    }
    else
    {
      return(L);
    }
  }
  else
  {
    //==================create anew ring with extra variables================

    execute("ring R1=("+charstr(bsr)+"),("+varstr(bsr)+",@y(1..m)),(dp(n),dp(m));");
    execute("minpoly=number("+mp+");");
    ideal id=imap(bsr,id);
    ideal A;

    for(z=1;z<=m;z++)
    {
      A[z]=lead(id[z])-@y(z);
    }

    A=groebner(A);
    ideal kern=nselect(A,1..n);// "kern" is the kernel of the ring map:
                        // R1----->bsr ;y(z)----> lead(id[z]).
                        //"kern" is the ideal of algebraic relations between
                        // lead(id[z]).
    export kern,A;// we export:
                  // * the ideal A  to avoid useless computations
                  //   between 2 steps in sagbi procedure.
                  // * the ideal kern : some times we can get intresting
                  //   informations from this ideal.
    setring bsr;
    map phi=R1,vars,id;

    // the sagbiSPolynomials are the image by phi of the generators of kern

    P=simplify(phi(kern),1);
    if(#[1]==0)
    {
      return(P);
    }
    else
    {
      L=P,R1;
      kill phi,vars;

      dbprint(printlevel-voice+3,"
// 'sagbiSPoly' created a ring as 2nd element of the list.
// The ring contains the ideal 'kern'  of algebraic relations between the
//leading terms of the generators of I.
// To access to this ring and see 'kern' you should give the ring a name,
// e.g.:
               def S = L[2]; setring S; kern;
      ");
      return(L);
    }
  }
}
example
{ "EXAMPLE:"; echo = 2;
 ring r=0, (x,y),dp;
 poly f1,f2,f3,f4=x2,y2,xy+y,2xy2;
 ideal I=f1,f2,f3,f4;
 sagbiSPoly(I);
 list L=sagbiSPoly(I,1);
 L[1];
 def S= L[2]; setring S; kern;
}

///////////////////////////////////////////////////////////////////////////////
static proc std1(ideal J,ideal I,list #)
  // I is contained in J, and it is assumed to be a standard bases!
  //This procedure computes a Standard basis for J from I one's
  //This procedure is essential for Spoly1 procedure.
{
  def br=basering;
  int tt;
  ideal Res,@result;

  if(size(#)>0) {tt=#[1];}

  if(size(I)==0) {@result=groebner(J);}

  if((size(I)!=0) && (size(J)-size(I)>=1))
  {
    qring Q=I;
    ideal J=fetch(br,J);
    J=groebner(J);
    setring br;
    Res=fetch(Q,J);// Res contains the generators that we add to I
                   // to get the generators of std(J);
    @result=Res+I;
  }

  if(tt==0) { return(@result);}
  else      { return(Res);}
}

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

static proc Spoly1(list l,ideal I,ideal J,int a)
  //an implementation of SAGBI construction Algorithm using Spoly
  //procedure leads to useless computations and affect the efficiency
  //of SAGBI bases computations. This procedure is a variant of Spoly
  //in order to avoid these useless compuations.
{
  degBound=0;
  def br=basering;
  ideal vars=maxideal(1);
  ideal B=ideal(br);
  int b=size(B);

  if(b!=0)
  {
    I=reduce(I,groebner(0));
    J=reduce(J,groebner(0));
  }
  int n,ii,jj=nvars(br),ncols(I),ncols(J);
  int z;
  list @L;
  string mp =string(minpoly);

  if(size(J)==0)
  {
    @L =sagbiSPoly(I,1);
  }
  else
  {
    ideal @sum=I+J;
    ideal P1;
    ideal P=l[1];//P is the ideal of spolynomials of I;
    def R=l[2];setring R;int kk=nvars(R);
    ideal J=fetch(br,J);

    //================create a new ring with extra variables==============
    execute("ring R1=("+charstr(R)+"),("+varstr(R)+",@y((ii+1)..(ii+jj))),(dp(n),dp(kk+jj-n));");
    // *levandov: would it not be easier and better to use
    // ring @Y = char(R),(@y((ii+1)..(ii+jj))),dp;
    // def R1 = R + @Y;
    // setring R1;
    // -> thus
    ideal kern1;
    ideal A=fetch(R,A);
    attrib(A,"isSB",1);
    ideal J=fetch(R,J);
    ideal kern=fetch(R,kern);
    ideal A1;
    for(z=1;z<=jj;z++)
    {
      A1[z]=lead(J[z])-var(z+kk);
    }
    A1=A+A1;
    ideal @Res=std1(A1,A,1);// the generators of @Res are whose we have to add
                            // to A to get std(A1).
    A=A+@Res;
    kern1=nselect(@Res,1..n);
    kern=kern+kern1;
    export kern,kern1,A;
    setring br;
    map phi=R1,vars,@sum;
    P1=simplify(phi(kern1),1);//P1 is th ideal we add to P to get the ideal
                              //of Spolynomials of @sum.
    P=P+P1;

    if (a==1)
    {
      @L=P,R1;
      kill phi,vars;
      dbprint(printlevel-voice+3,"
// 'Spoly1' created a ring as 2nd element of the list.
// The ring contains the ideal 'kern'  of algebraic relations between the
//generators of I+J.
// To access to this ring and see 'kern' you should give the ring a name,
// e.g.:
               def @ring = L[2]; setring @ring ; kern;
      ");
    }
    if(a==2)
    {
      @L=P1,R1;
      kill phi,vars;
    }
  }
  return(@L);
}
///////////////////////////////////////////////////////////////////////////////

proc sagbiReduction(poly p,ideal dom,list #)
"USAGE: sagbiReduction(p,dom[,n]); p poly , dom  ideal
RETURN: a polynomial, after one step subalgebra reduction
@format
    Three algorithm variants are used to perform subalgebra reduction.
    The positive interger n determines which variant should be used.
    n may take the values 0 (default), 1 or 2.
@end format
EXAMPLE: sagbiReduction; shows an example"
{
  def bsr=basering;
  ideal B=ideal(bsr);//When the basering is quotient ring  this type casting
                     // gives the quotient ideal.
  int b=size(B);
  int n=nvars(bsr);

  //In quotient rings, SINGULAR, usually does not reduce polynomials w.r.t the
  //quotient ideal,therefore we should first  reduce ,when it is necessary for computations,
  // to have a uniquely determined representant for each equivalent
  //class,which is the case of this algorithm.

  if(b !=0) //means that the basering is a quotient ring
  {
    p=reduce(p,std(0));
    dom=reduce(dom,std(0));
  }

  int i,choose;
  int z=ncols(dom);

  if((size(#)>0) && (typeof(#[1])=="int"))
  {
    choose = #[1];
  }
  if (size(#)>1)
  {
    choose =#[2];
  }

  //=======================first algorithm(default)=========================
  if ( choose == 0 )
  {
    list L = algebra_containment(lead(p),lead(dom),1);
    if( L[1]==1 )
    {
      // the ring L[2] = char(bsr),(x(1..nvars(bsr)),y(1..z)),(dp(n),dp(m)),
      // contains poly check s.t. LT(p) is of the form check(LT(f1),...,LT(fr))
      def s1 = L[2];
      map psi = s1,maxideal(1),dom;
      poly re = p - psi(check);
      // divide by the maximal power of #[1]
      if ( (size(#)>0) && (typeof(#[1])=="poly") )
      {
        while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
        {
          re=re/#[1];
        }
      }
      return(re);
    }
    return(p);
  }
  //======================2end variant of algorithm=========================
  //It uses two different commands for elimaination.
  //if(choose==1):"elimainate"command.
  //if (choose==2):"nselect" command.
  else
  {
    poly v=product(maxideal(1));

    //------------- change the basering bsr to bsr[@(0),...,@(z)] ----------
    execute("ring s=("+charstr(basering)+"),("+varstr(basering)+",@(0..z)),dp;");
    // Ev hier die Reihenfolge der Vars aendern. Dazu muss unten aber entsprechend
    // geaendert werden:
    //  execute("ring s="+charstr(basering)+",(@(0..z),"+varstr(basering)+"),dp;");

    //constructs the leading ideal of dom=(p-@(0),dom[1]-@(1),...,dom[z]-@(z))
    ideal dom=imap(bsr,dom);
    for (i=1;i<=z;i++)
    {
      dom[i]=lead(dom[i])-var(nvars(bsr)+i+1);
    }
    dom=lead(imap(bsr,p))-@(0),dom;

    //---------- eliminate the variables of the basering bsr --------------
    //i.e. computes dom intersected with K[@(0),...,@(z)].

    if(choose==1)
    {
      ideal kern=eliminate(dom,imap(bsr,v));//eliminate does not need a
                                            //standard basis as input.
    }
    if(choose==2)
    {
      ideal kern= nselect(groebner(dom),1..n);//"nselect" is combinatorial command
                                         //which uses the internal command
                                         // "simplify"
    }

    //---------  test wether @(0)-h(@(1),...,@(z)) is in ker ---------------
    // for some poly h and divide by maximal power of q=#[1]
    poly h;
    z=size(kern);
    for (i=1;i<=z;i++)
    {
      h=kern[i]/@(0);
      if (deg(h)==0)
      {
        h=(1/h)*kern[i];
        // define the map psi : s ---> bsr defined by @(i) ---> p,dom[i]
        setring bsr;
        map psi=s,maxideal(1),p,dom;
        poly re=psi(h);
        // divide by the maximal power of #[1]
        if ((size(#)>0) && (typeof(#[1])== "poly") )
        {
          while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
          {
            re=re/#[1];
          }
        }
        return(re);
      }
    }
    setring bsr;
    return(p);
  }
}
example
{"EXAMPLE:"; echo = 2;
  ring r= 0,(x,y),dp;
  ideal dom =x2,y2,xy-y;
  poly p=x4+x3y+xy2-y2;
  sagbiReduction(p,dom);
  sagbiReduction(p,dom,1);
  sagbiReduction(p,dom,2);
}

///////////////////////////////////////////////////////////////////////////////
static proc completeReduction(poly p,ideal dom,list#)//reduction
{
  poly p1=p;
  poly p2=sagbiReduction(p,dom,#);
  while (p1!=p2)
  {
    p1=p2;
    p2=sagbiReduction(p1,dom,#);
  }
  return(p2);
}
///////////////////////////////////////////////////////////////////////////////

static proc completeReduction1(poly p,ideal dom,list #) //tail reduction
{
  poly p1,p2,re;
  p1=p;
  while(p1!=0)
  {
    p2=sagbiReduction(p1,dom,#);
    if(p2!=p1)
    {
      p1=p2;
    }
    else
    {
      re=re+lead(p2);
      p1=p2-lead(p2);
    }
  }
  return(re);
}



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

proc sagbiNF(id,ideal dom,int k,list#)
"USAGE: sagbiNF(id,dom,k[,n]); id either poly or ideal,dom ideal, k and n positive intergers.
RETURN: same as type of id; ideal or polynomial.
@format
    The integer k determines what kind of s-reduction is performed:
    - if (k=0) no tail s-reduction is performed.
    - if (k=1) tail s-reduction is performed.
    Three Algorithm variants are used to perform subalgebra reduction.
    The positive integer n determines which variant should be used.
    n may take the values (0 or default),1 or 2.
@end format
NOTE: computation of subalgebra normal forms may be performed in polynomial rings or quotients
      thereof
EXAMPLE: example sagbiNF; show example "
{
  int z;
  ideal Red;
  poly re;
  if(typeof(id)=="ideal")
  {
    int i=ncols(id);
    for(z=1;z<=i;z++)
    {
      if(k==0)
      {
        id[z]=completeReduction(id[z],dom,#);
      }
      else
      {
        id[z]=completeReduction1(id[z],dom,#);//tail reduction.
      }
    }
    Red=simplify(id,7);
    return(Red);
  }
  if(typeof(id)=="poly")
  {
    if(k==0)
    {
      re=completeReduction(id,dom,#);
    }
    else
    {
      re=completeReduction1(id,dom,#);
    }
    return(re);
  }
}
example
{"EXAMPLE:"; echo = 2;
 ring r=0,(x,y),dp;
 ideal I= x2-xy;
 qring Q=std(I);
 ideal dom =x2,x2y+y,x3y2;
 poly p=x4+x2y+y;
 sagbiNF(p,dom,0);
 sagbiNF(p,dom,1);// tail subalgebra reduction is performed
}


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

static proc intRed(id,int k, list #)
{
  int i,z;
  ideal Rest,intRed;
  z=ncols(id);
  for(i=1;i<=z;i++)
  {
    Rest=id;
    Rest[i]=0;
    Rest=simplify(Rest,2);
    if(k==0)
    {
      intRed[i]=completeReduction(id[i],Rest,#);
    }
    else
    {
      intRed[i]=completeReduction1(id[i],Rest,#);
    }
  }
  intRed=simplify(intRed,7);//1+2+4 in simplify command
  return(intRed);
}
//////////////////////////////////////////////////////////////////////////////

proc sagbi(id,int k,list#)
"USAGE: sagbi(id,k[,n]); id ideal, k and n positive integers.
RETURN: A SAGBI basis for the subalgebra defined by the generators of id.
@format
    k determines what kind of s-reduction is performed:
     - if (k=0) no tail s-reduction is performed.
     - if (k=1) tail s-reduction is performed, and S-interreduced SAGBI basis
       is returned.
    Three algorithm variants are used to perform subalgebra reduction.
    The positive interger n determine which variant should be used.
    n may take the values (0 or default),1 or 2.
@end format
NOTE: SAGBI bases computations may be performed in polynomial rings or quotients
      thereof.
EXAMPLE: example sagbi; show example "
{
  degBound=0;
  ideal S,oldS,Red;
  list L;
  S=intRed(id,k,#);
  while(size(S)!=size(oldS))
  {
    L=Spoly1(L,S,Red,2);
    Red=L[1];
    Red=sagbiNF(Red,S,k,#);
    oldS=S;
    S=S+Red;
  }
  return(S);
}
example
{ "EXAMPLE:"; echo = 2;
 ring r= 0,(x,y),dp;
 ideal I=x2,y2,xy+y;
 sagbi(I,1,1);
}
///////////////////////////////////////////////////////////////////////////////
proc sagbiPart(id,int k,int c,list #)
"USAGE: sagbiPart(id,k,c[,n]); id ideal, k, c and n positive integers
RETURN: A partial SAGBI basis for the subalgebra defined by the generators of id.
@format
     k determines what kind of s-reduction is performed:
     - if (k=0) no tail s-reduction is performed.
     - if (k=1) tail s-reduction is performed, and S-intereduced SAGBI basis
      is returned.
     c determines, after how many loops the Sagbi basis computation should stop.
     Three algorithm variants are used to perform subalgebra reduction.
     The positive integer n determines which variant should be used.
     n may take the values (0 or default),1 or 2.
@end format
NOTE:- SAGBI bases computations may be performed in polynomial rings or quotients
       thereof.
     - This version of sagbi is interesting in the case of subalgebras with infinte
       SAGBI basis. In this case, it may be used to check, if the elements of this
       basis have a particular form.
EXAMPLE: example sagbiPart; show example "
{
  degBound=0;
  ideal S,oldS,Red;
  int counter;
  list L;
  S=intRed(id,k,#);
  while((size(S)!=size(oldS))&&(counter<=c))
  {
    L=Spoly1(L,S,Red,2);
    Red=L[1];
    Red=sagbiNF(Red,S,k,#);
    oldS=S;
    S=S+Red;
    counter=counter+1;
  }
  return(S);
}
example
{ "EXAMPLE:"; echo = 2;
 ring r= 0,(x,y),dp;
 ideal I=x,xy-y2,xy2;//the corresponding Subalgebra has an infinte SAGBI basis
 sagbiPart(I,1,3);// computations should stop after 3 turns.
}
//////////////////////////////////////////////////////////////////////////////