LeadingMonomial[poly_,vars_]:=
Module[{p,exps},
p=Expand[poly];
If[poly===0,Throw["The zero polynomial does not have a leading monomial."]];
If[Head[p]===Plus,p=List@@p,p={p}];
exps=Exponent[#,vars]&/@p;
exps=Sort[exps];
Return[Times@@(vars^Last[exps])];
];

LeadingCoefficient[poly_,vars_]:=
Module[{lm,c},
lm=LeadingMonomial[poly,vars];
c=poly;
Do[
c=Coefficient[c,vars[[i]],Exponent[lm,vars[[i]]]];
,{i,1,Length[vars]}];
Return[c];
];

PolynomialReduction[poly_,set_,vars_]:=
Module[{lm,lc,rem=0,p=poly,i,test},
While[p=!=0,
lm=LeadingMonomial[p,vars];
lc=LeadingCoefficient[p,vars];
i=1;
While[i≤Length[set],
test=lm/LeadingMonomial[set[[i]],vars];
If[PolynomialQ[test,vars],
p=Expand[p-lc/LeadingCoefficient[set[[i]],vars]*test*set[[i]]];
Break[];
];
i++;
];
If[i>Length[set],rem=rem+lc*lm;p=p-lc*lm];
];
Return[rem];
];

p1=x*y-1;
p2=x^2+y+1;

PolynomialReduction[x^2*y,{p1,p2},{x,y}]

PolynomialReduction[x^2*y,{p2,p1},{x,y}]

(* Buchberger's algorithm computes a Groebner bases for the ideal *)
Buchberger[set1_,vars_]:=
Module[{i,pairs,set=set1,p1,p2,lm1,lm2,gcd,sp},
pairs=Subsets[set,{2}];
While[pairs=!={},
p1=pairs[[1,1]];
p2=pairs[[1,2]];
pairs=Rest[pairs];
lm1=LeadingMonomial[p1,vars];
lm2=LeadingMonomial[p2,vars];
gcd=PolynomialGCD[lm1,lm2];
sp=Expand[lm2/gcd*LeadingCoefficient[p2,vars]*p1-lm1/gcd*LeadingCoefficient[p1,vars]*p2];
sp=PolynomialReduction[sp,set,vars];
If[sp=!=0,
pairs=Join[pairs,Table[{set[[i]],sp},{i,Length[set]}]];
AppendTo[set,sp];
];
];
Return[set];
];

Examples

Find common zeros of a polynomial system

Gröbner bases can be used for solving a polynomial system. When using a lexicographic monomial order, an elimination is performed. Hence we get a univariate equation which we can solve, and backwards substitute the results in the other equations:

Buchberger[{x*y-1,x^2+y+1},{x,y}]

Solve[Last[%]==0,y]

Solve[#==0,x]&/@((x*y-1)/.%)

Uniqueness of remainder

ideal={x*y-1,x^2+y+1}

PolynomialReduction[x^2*y,ideal,{x,y}]

PolynomialReduction[x^2*y,Reverse[ideal],{x,y}]

gb=Buchberger[ideal,{x,y}]

PolynomialReduction[x^2*y,gb,{x,y}]

PolynomialReduction[x^2*y,Reverse[gb],{x,y}]

gcd in the univariate case

Expand[(x-3)*(x+1)]

Buchberger[{x^3-10x+3,x^2-2x-3},{x}]

PolynomialGCD[x^3-10x+3,x^2-2x-3]

Gauss in the first-degree case

Buchberger[{a+b+c-1,a-2b-3c+4,2a-b+2c-2},{a,b,c}]

Solve[{a+b+c-1==0,a-2b-3c+4==0,2a-b+2c-2==0},{a,b,c}]

Compare also with the Mathematica command GroebnerBasis!