/////// PROCEDURE TO COMPUTE SASBI-STANDARD BASES //////

// The following libraries will be used: 

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


////////// Procedure to comupte Sasbi-Standard Weak normal form ////////

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

//The following procedure needs two more procedures SSDL (provides a list of vectors "ag" such that lm(f)=lm(ag) with a in A and g in I) 
and SortE (provides a sorted ascending list of vectors from SSDL with respect to ecart). 

// Main Procedure:
proc ModWSSNF(vector f, module I, ideal A)

"USAGE: ModWSSNF(f,I,A); f a polynomial vector, I an A-module, A a subalgebra (which has a finite Sasbi basis).

RETURN: a polynomial vector h."

{
   A=Sasbi(A);
   module G=I;
   vector h=f;
   vector h1;
   poly c;
   list D;
   map psi ;
   int i,k;
   while(h!=0 && h1!=h)
    {  
       D=SSDL(h,G,A);
       D=sortE(D);
       if (size(D)>0)	
          {		
	    vector hi=D[1];
          if (ecart(h)<ecart(hi))   
            {
             G[size(G)+1]=h;
            }   
	    c=leadcoef(h)/leadcoef(hi);
          h1=h;
          h=h-c*hi; 
	    kill hi;
	    }
       if(size(D)==0)
         {
           h1=h;  
         }
     }
   return(h);
} 

// Procedure 1:

proc SSDL(vector f, module I, ideal A)

"USAGE: SSDL(f,I,A); f a polynomial vector, I an A-module, A a subalgebra (which has a finite Sasbi basis).

 RETURN: a list M."

{
   A=sagbi(A);
   int u;
   u=nrows(I);
   module G=I;
   vector h=f ;
   poly q,j;
   list L,M ;
   int i,k,t;
   t=1;
   vector h1=lead(h);
   for(k=1;k<=size(G);k++)
    {
     vector gk=lead(G[k]);
     poly sumg,sum1;
     for(i=1;i<=u;i++)
     {
     sumg=sumg+gk[i];
     sum1=sum1+h1[i];
     }
     q=sum1/sumg;
     vector p=q*gk;
     kill gk,sumg,sum1;
     if(q!=0 && h1==p)
     {
     map psi;
     L= algebra_containment(q,lead(A),1); 
     if(L[1]==1)
     {
      def s=L[2];
      psi= s,maxideal(1),A ;
      j=psi(check);
      M[t]=j*G[k];
      kill s;
      t=t+1;
      }
     }
    kill p;
    }
  return(M);
}

// Procedure 2:

 proc sortE(list L)

"USAGE: sortE(L); L a list of vectors.

 RETURN: a list M."

{ 
 int i,j;
 int n=size(L);
 list M=L;
 for(j=1;j<=n;j++)
    { 
     for(i=j+1;i<=n;i++)
        { 
          int hj=ecart(M[j]);
          int hi=ecart(M[i]);
          if (hi<hj)
	    {
	      vector Hi=M[i];
	      vector Hj=M[j];
            M[j]=Hi;
            M[i]=Hj;
           }
          //kill hi,hj,Hi,Hj;
        } 
    }
 return(M);
}


//// Procedure to compute a generating set for syzygies.////

// For a subalgebra A with a finite Sasbi basis and an A-module I, the following procedure computes the 
leading term generating set for syzygies of lead(I) over lead(A). It needs a procedure ModSYZ for computing syzygies through Grobner basis technique.

//Main Procedure:

proc ModLTSS(module I,ideal A)

"USAGE: ModLTSS(I,A); I an A-module, A a subalgebra (which has a finite Sasbi basis).

RETURN: a module M."

{
  def r=basering;
  A=Sasbi(A);
  int n=nvars(r);
  int m=size(A);
  int k=size(I);
  int p=nrows(I);
  ideal B=lead(A);
  module J=lead(I);
  ideal vars = maxideal(1);
  int c;
  module Q1 ;
  if(I==0)
  {
   return(Q1);
  }
  else
  {
   execute("ring R1=("+charstr(r)+"),(y(1..m),"+varstr(r)+"),(ds(m),ds(n));");
   list G;
   ideal B=imap(r,B);
   module J=imap(r,J);
   ideal A=imap(r,A);
   int i,j;
   module M;
   for(i=1;i<=k;i++)
    { 
     poly sumgi;
     vector Ji=J[i];
     for(j=1;j<=p;j++)
     {
     sumgi=sumgi+J[i][j];
     }
     list Li=algebra_containment(sumgi,B,1);
     if (Li[1]==1)
     {	
       def Si=Li[2];
       poly ji=imap(Si,check);
       G[i]=ji;
     }
vector qi=J[i]/sumgi;
vector vi;
for(c=1;c<=size(G);c++)
{
poly fc=G[c];
vi=fc*qi; 
}
M[i]=vi;
} //Lead of I in PP of lead of A
 module H=std(M);
 list N=ModSYZ(M);
 execute("ring R=("+charstr(r)+"),(y(1..m)),(ds(m));");
 setring r;
 map phi= R,B;
 setring R;
 ideal T=kernel(r,phi);
 setring R1;
 ideal T=imap(R,T);
 int i;
 module KER;
 for(i=1;i<=p;i++)
 {
  KER=KER+T*gen(i);
 }
 module J=intersect(H,KER);
 list Q;
 int b=size(J);
 int c=size(H);
 if(b!=0)
 {
  int i,d,e;
  list D;
  for(d=1;d<=b;d++)
  {
   vector qd; 
   vector gd=J[d]; 
   list td=division(gd,H);
   int sd;
   for(e=1;e<=c;e++)
   {
   for(i=1;i<=c;i++)
   {
   if(size(td[1][i,1])==0)
    {
     sd=sd+1;
    }
   else
    {
     sd=sd;
    }
   }
   int s=k-sd;
   if(s!=1) 
    {
     vector pde=td[1][e,1]*gen(e);
     qd=qd+pde;
    }
    else
    {
     qd=qd;     
     }
    }
     D[d]=qd;
     }
     Q=N,D; 
     }
     else
     {
     Q=N,0;
     }
     if(size(Q[1])==0 && size(Q[2])==0)
     {
     setring r;
     module M;
     return(M);
     }
     else
     {
     setring r;
     map psi=R1,A,maxideal(1);
     list W=psi(Q);
     int i,j,l,t;
     module M;
     t=1;
     int n=size(W[t]);
     int f=1;
     module Z;
     for(i=1;i<=n;i++)
     {
     for(j=i+1;j<=size(W[t][i]);j++)
      {
       vector vij=W[t][i][j];
       Z[f]=vij;
       f=f+1;
      }
      M=Z;
     }
      t=t+1;
	for(l=1;l<=size(W[2]);l++)
	{
       vector vt=W[t][l];
       if(vt!=0)
        {
         M=M,vt;
         }
       } 
      return(M);
      }
     }
 }

// Procedure 1:

proc ModSYZ(module G)

"USAGE: ModSYZ(G), G a module.

RETURN: a list T."

{
  int i,j,l,k;
  int u=nrows(G);
  int l=size(G);
  list T;
  for(i=1;i<=l-1;i++)
    { 
     poly sumgi;
     vector gi=G[i];
     for(k=1;k<=u;k++)
     {
     sumgi=sumgi+gi[k];
     }
     vector qi=gi/sumgi;
     list S;
     for(j=i+1;j<=l;j++)
	 {
      poly sumgj;
      vector gj=G[j];
      for(k=1;k<=u;k++)
       {
       sumgj=sumgj+gj[k];
       }
       vector qj=gj/sumgj;
	 if(qi==qj)
          {
           poly gij=lcm(sumgi,sumgj);
          }
       else
          {
           poly gij=0;
          } 
	     vector vij=(gij/sumgi)*gen(i)-(gij/sumgj)*gen(j);
	     S[j]=vij;
        } 
        T[i]=S;
      }
  return(T);
}


//// Procedure to compute Sasbi-Standard basis //// 

//The following procedure computes a Sasbi-Standard basis of an A-module I over sublagebra A. 
It needs two procedures ModSPOLY (computes S polynomial vectors which further uses the procedure ModLTSS) and ModWSSNF (defined earlier).

// Main Procedure:

proc ModSSB(module I,ideal A)

"USAGE: ModSSB(I,A); I an A-module, A a subalgebra (which has a finite Sasbi basis).

 RETURN: a module MSSB."

{
  module MSSB,oldMSSB;
  vector Red;
  module T;
  int k,l;
  MSSB=I;
  while(size(MSSB)!=size(oldMSSB))
    {
	A=Sasbi(A);
	T=ModSPOLY(MSSB,A);      
	l=size(T);
	if(l==0)
	{
	 oldMSSB=MSSB;
	}
	else
	{
	for(k=1;k<=l;k++)
       {
	 Red=T[k];
	 Red=ModWSSNF(Red,MSSB,A);
	 oldMSSB=MSSB;
	 MSSB=MSSB+Red;
	 }
      }
    }
  return(MSSB) ;
}

// Prodcedure 1:

proc ModSPOLY(module I,ideal A)

"USAGE: ModSPOLY(I,A); I an A-module, A a subalgebra.

 RETURN" a module S."

{
 int i,j;
 module P=ModLTSS(I,A);
 module S;
 for(j=1;j<=size(P);j++)
 {
  vector Vj=P[j];
  vector hj;
  for(i=1;i<=size(I);i++)
   {
    vector pij = Vj[i]*I[i];
    hj=hj+pij;
   }
    S=S+hj;
  }
  return(S);
}

/// Singular Examples for all main procedures:

//Example 1:

ring r=0,(z,y,x),ls;
ideal A=x2 + y, y4, z;
module I=[x2 + y,y4],[y8 + 2y4, y6];
vector f=[x2 + z, z3];
ModWSSNF(f,I,A);
-y*gen(1)-y4*gen(2)+z*gen(1)+z3*gen(2)
ModLTSS(I,A);
_[1]=-x2*gen(2)-y*gen(2)+2y4*gen(1)
ModSPOLY(I,A);
_[1]=-y6x2*gen(2)-y7*gen(2)+2y8*gen(2)-y8x2*gen(1)-y9*gen(1)
ModSSB(I,A);
_[1]=x2*gen(1)+y*gen(1)+y4*gen(2)
_[2]=2y4*gen(1)+y6*gen(2)+y8*gen(1)
_[3]=-y6x2*gen(2)-y7*gen(2)+2y8*gen(2)-y8x2*gen(1)-y9*gen(1)

//Example 2:

ring r=0,(x,y),Ds;
module I=[x,y2],[y2,y3];
ideal A=x,y2,y3;  
vector f=[x2,x4];
ModWSSNF(f,I,A);
-xy2*gen(2)+x4*gen(2)
ModLTSS(I,A);
_[1]=-x*gen(2)+y2*gen(1)
ModSPOLY(I,A);
_[1]=-xy3*gen(2)+y4*gen(2)
ModSSB(I,A);
_[1]=x*gen(1)+y2*gen(2)
_[2]=y2*gen(1)+y3*gen(2)
_[3]=-xy3*gen(2)+y4*gen(2)

//Example 3:

ring r=0,(z,y,x),ds;
module I=[z2+2zy+y2,z+y],[xy2,xy2z+xy3];
ideal A=z+y,xy2;
vector f=[x2y4,z+y];
ModWSSNF(f,I,A);
-z2*gen(1)-2zy*gen(1)-y2*gen(1)+y4x2*gen(1)
ModLTSS(I,A);
_[1]=0
ModSPOLY(I,A);
_[1]=0
ModSSB(I,A);
_[1]=z*gen(2)+y*gen(2)+z2*gen(1)+2zy*gen(1)+y2*gen(1)
_[2]=y2x*gen(1)+zy2x*gen(2)+y3x*gen(2)

//Example 4:

ring r=0,(x,y),Ds; 
module I=[x,y2],[y2,xy];
ideal A=x,y2,xy;
ModLTSS(I,A);
_[1]=0
ModSPOLY(I,A);
_[1]=0
ModSSB(I,A);
_[1]=x*gen(1)+y2*gen(2)
_[2]=xy*gen(2)+y2*gen(1)