TermePrincipal := function (P)
  return Coefficients(P)[1]*Monomials(P)[1];
end function;

ChercheMultiple := function(M, TPI)
  local p, i;
  for i:=1 to #TPI do
    if IsDivisibleBy(M,TPI[i]) then
       return true,i;
    end if;
  end for;
  return false,_;
end function;

Division := function(AP, I)
  local P, TPI, L, R, M, t, i, q;
  P   := AP;
  TPI := [TermePrincipal(i) : i in I];
  L   := [0*x : x in I];
  R   := 0*P;
  while not IsZero(P) do
    M := TermePrincipal(P);
    t, i := ChercheMultiple(M, TPI);
    if t then
      q     := M div TPI[i];
      L[i] +:= q;
      P    -:= q * I[i];
    else 
      R +:= M;
      P -:= M;
    end if;
  end while;
  return L, R;
end function;

Syzygy:= function(gi,gj)
  local d;
  d:=GCD(gi,gj);
  return gi div d, gj div d;
end function;

CheckGrob := function(B)
  local TPB, i, j, mij, mji, L ,h;
  TPB := [TermePrincipal(b) : b in B];
  for i:=1 to #B do
    for j:=i+1 to #B do
      mij, mji := Syzygy(TPB[i],TPB[j]);
      L, h := Division(mji*B[i]-mij*B[j],B);
      if not IsZero(h) then
        return false, h;
      end if;
    end for;
  end for;
  return true,_;
end function;

Grob := function(AB)
  local B, t, h;
  B:=AB;
  repeat
    t, h := CheckGrob(B);
    if t then
      return B;
    end if;
    Append(~B, h);
  until false;
end function;

Nullstellensatz:= function(P, B)
  local e, L, R, Q;
  e:=1; Q:=P;
  repeat 
    L, R := Division (Q, B);
    if IsZero(R) then
      return e;
    end if;
    e +:= 1; Q *:= P;
  until false;
end function;

A<X11,X12,X21,X22> := PolynomialAlgebra(RationalField(),4);

M := Matrix(A, 2, 2, [X11,X12,X21,X22]);
M2 := M^2;
I := [M2[1,1],M2[1,2],M2[2,1],M2[2,2]];
B:=Grob(I);
Nullstellensatz(Trace(M),B);
Nullstellensatz(Determinant(M),B);

A<X,Y> := PolynomialAlgebra(RationalField(),2);
B:=Grob([X^2,X*Y+Y^2]);
Nullstellensatz(X+Y,B);

A<X,Y,Z> := PolynomialAlgebra(RationalField(),3);
Grob([Y^2,X*Y*Z-Z^3]);

