/////// PROCEDURE TO COMPUTE SASBI BASES //////

// The following libraries will be used: 

LIB"algebra.lib";
LIB"sagbi.lib";
LIB "grobcov.lib";
LIB"elim.lib" ; //for "nselect" procedure
LIB"teachstd.lib"; //for ecart

////////// Procedure to compute Sasbi Weak normal form ////////

// Let A=B_> be a localization of a polynomial subalgebra B with respect to a local monomial ordering >. 
For a polynomial f of R_> (a localization of ring R with respect to >) and a finite set of polynomials I in A, 
the following procedure computes a weak Sasbi normal form of f with respect to A.

proc WSNF(poly f,ideal I)

"USAGE: WSNF(f,I); f a polynomial, I a subalgebra.

RETURN: a polynomial h."

{ 
  ideal G=I;
  poly h=f;
  poly h1,j;
  list L;
  map psi;
  while(h!=0 && h1!=h) 
  {
     L= algebra_containment(lead(h),lead(G),1);
     if (L[1]==1)
     {
        def s= L[2];
        psi= s,maxideal(1),G; 
        j= psi(check);
        if (ecart(h)<ecart(j))   
        {
           G[size(G)+1]=h;
        }          
        h1=h;          
        h=h-j;
        kill s;
     }    
  if (L[1]==0){h1=h;}
   }
   return (h);
 }


//// Procedure to compute Sasbi-Spolynomials ////

// For a subalgebra A, the following procedure computes S polynomials of A.

proc sasbiSpoly(ideal id)

"USAGE: sasbiSpoly(id); id a subalgebra.

 RETURN: an ideal p."

{
 def bsr= basering ;
 ideal vars = maxideal(1) ;
 int n=nvars(bsr) ;
 int m=ncols(id) ;
 int z ;
 ideal p ;
 if(id==0)
 {
  return(p) ;
 }
 else
 {
  execute("ring R1=("+charstr(bsr)+"),(@y(1..m),"+varstr(bsr)+"),(ds(m),ds(n));");
  ideal id =imap(bsr,id) ;
  ideal A ;
  for (z=1; z<=m; z++)
  {
   A[z]=lead(id[z])-@y(z) ;
  }
  A=std(A) ;
  ideal kern=nselect(A,m+1..m+n) ;
  export kern,A ;
  setring bsr ;
  map phi= R1,id ;
  p=simplify(phi(kern),2) ;
  return (p) ;
 }
}

//// Procedure to compute Sasbi basis //// 

//The following procedure computes a Sasbi basis of a sublagebra A. 
It needs two procedures sasbiSpoly and WSNF (both are defined earlier).

proc Sasbi(ideal id)

"USAGE: Sasbi(id); id a subalgebra.

 RETURN: an ideal S."

{
 ideal S,oldS,Red ;
 list L ;
 int z,n ;
 S=id ;
 while( size(S)!=size(oldS))
 {
  L=sasbiSpoly(S) ;
  if (L[1][1]==0) 
  {return (S);}
  n=size(L[1]) ;
  for (z=1; z<=n; z++)
  {
   Red=L[1][z] ;
   Red=WSNF(Red[1],S) ;
   oldS=S ;
   S=S+Red ;
  }
 }
return(S) ;
}


/// Singular Examples for all main procedures:

//Example 1:
ring r=0,(z,y,x),ls;
ideal A=x2 + y, y4, z;
poly f=x2+z;
WSNF(f,A);
-y+z
sasbiSpoly(A);
_[1]=0
Sasbi(A);
_[1]=x2+y
_[2]=y4
_[3]=z

//Example 2:
ring r=0, (x,y), ls;
ideal A=x2+y8, x4+x6, x2y4+x4;
poly f=y16+x2;
WSNF(f,A);
x2-2x2y8-x4
sasbiSpoly(A);
_[1]=-y16-2x2y8+x6
_[2]=-y16+x2y4-2x2y8
Sasbi(A);
_[1]=y8+x2
_[2]=x4+x6
_[3]=x2y4+x4

//Example 3:
ring r=0, (x,y,z), Ds;
ideal A= x2+y, y2+x6, z2;
poly f= -20y2+x4;
WSNF(f,A);
x4+20x6
sasbiSpoly(A);
_[1]=-2x2y-x4+x6
Sasbi(A);
_[1]=y+x2
_[2]=y2+x6
_[3]=z2
_[4]=-2x2y-x4+x6

//Example 4
ring r=0, (x), ds;
ideal A=x8,x10+x13,x12+x15;
poly f= x20+x28;
WSNF(f,A);
-x23+x28
sasbiSpoly(A);
_[1]=x23+x26
_[2]=2x27+x30
Sasbi(A);
_[1]=x8
_[2]=x10+x13
_[3]=x12+x15
_[4]=x23+x26
_[5]=2x27+x30