-- We define three functions, one of which - Buchberger - shall be 
-- used to compute a G.b. for the ideal generated by polynomials F 

Buchberger = F -> (
     local G,S,i,j,p,f; -- make all the variables used in the routine local 
     -- initialize the current basis
     G = F; 
     -- initialize the queue of s-pairs
     S = {};
     for i from 0 to #G-1 do
     for j from i+1 to #G-1 do
          S = S  | {(G#i,G#j)};
     -- main loop
     while #S>0 do (
	  p = first S; S = drop(S,1); -- pick the first s-pair in S
	  f = reduce(Spoly(p#0,p#1), G);  
	  if f != 0 then (
	       -- "apply" applies a function (second argument)
	       -- to a list (first)
	       S = S | apply(G, g->(g,f));
	       G = G | {f} 
	       )     
	  );
     G
     )

reduce = (f,G) -> (
     ret := 0;
     while f!=0 do (
     	 while f!=0 and any(G, g->(leadTerm f)%(leadTerm g)==0) do (
	     -- select one element of G s.t. its LT divides the LT(f)
     	     g := select(1, G, g->(leadTerm f)%(leadTerm g)==0);
	     -- the result of the previous operation is a list of one element
	     f = f%(first g);   
	     );
	 if f!=0 then (
	     ret = ret + leadTerm f;
	     f = f - leadTerm f;
	     )
	 );
     ret
     )

Spoly = (f,g) -> (
     m := lcm(leadMonomial f, leadMonomial g);
     (m//leadTerm f)*f - (m//leadTerm g)*g
     )
end

restart
load "buchberger.m2"
R = QQ[x,y];
F = {x^3,x^2*y-y^3};
Buchberger F
gens gb ideal F -- compare with the built-in command "gb"

R = QQ[x,y,Weights=>{1,10}];
F = {x^3,x^2*y-y^3};
Buchberger F
gens gb ideal F -- compare with the built-in command "gb"

R = QQ[x,y,MonomialOrder=>Eliminate 1]; -- use elimination order w.r.t. x
F = {x^3,x^2*y-y^3};
Buchberger F
gens gb ideal F -- compare with the built-in command "gb"