# Now we revisit the process of  Grobner bases and finding solutions of
# a system of polynomial equations using MAPLE.
> restart; F:=[x^2+y^2-1,x^2-y^2-z,x+y-z^2];
> with(Groebner);
> G:=gbasis(F,plex(z,y,x));
> factor(G);
> _EnvExplicit:=true; solve(4*x^4+4*x^3-2*x^2-4*x-1,x);
# This environment variable will try to express with radicals all those
# roots for which is possible (e.g. in polynomials of degree less than
# or equal to four). When the variable is set to false then roots are
# handled symbolically.

> assign({solve(G[3],{z}), solve(G[2],{y})});
> rootlist:=[solve(G[1])];

> for x in rootlist do
>  whatever:=simplify(['x'=x,'y'=y, 'z'=z])):
>  if simplify(subs(op(whatever),F))=[0,0,0] then 
>    
> print("valid solution"); end if;
> od:
# We can also (A) decide whether a system has solutions
# (B) decide whether there are finitely many or infinitely many
# solutions
# (C) decide exactly how many solutions are there.

# LEMMA: (Hilbert's Nullstellensatz) 
# A system of polynomial equations has complex solution if 
# and only if a Grobner basis for the ideal generated by the 
# polynomials is not equal [1].

> restart; with(Groebner);
> gbasis([x+x*y^2-1,x^2*y+y-1,x^2+y^2-1],tdeg(x,y));
# LEMMA: A system of equations has a finite number of solutions if and
# only if any
# Grobner basis of the polynomials has the property that for every
# variable
# x_i there exists a polynomial such that its leading term with respect
# to the
# chosen term is a power of x_i.
> polys:=[c*x+x*y^2+x*z^2-1,c*y+y*x^2+y*z^2-1,c*z+z*x^2+z*y^2-1];
> with(Groebner): is_finite(gbasis(polys,tdeg(x,y,z)));
> sys:=gsolve(polys,[c,x,y,z]);

> for s in sys do solve({op(s[1])},{x,y,z}) end do;
> G:=gbasis(polys,plex(c,x,y,z));

> is_finite(G);



# What MAPLE uses to determine whether the number of solutions is finite
# is the following lemma:

# LEMMA: For a system of polynomials the number of solutions is finite
# if
# and only if the set of monomials which are not multiples of a the
# leading
# monomials of a Grobner basis is finite.

# Now, if the number of solutions is infinite there is a nice way we
# can know
# the dimension of the solution space (i.e. how man free parameters
# should
# there be?).


> hilbertdim(polys,tdeg(c,x,y,z));


# Suppose now we want to actually count the number of solutions:

# LEMMA: Suppose that the polynomials f_1, f_2, ..., f_s have a finite
# number of solutions. Then, the number of solutions (counted with
# multiplicities and solution at infinity) is equal to the cardinality
# of the sets of monomials that are no multiples of leading monomials
# of any Grobner basis (any term order).

> G:=gbasis(polys, tdeg(x,y,z));
> LM:=map(f -> leadmon(f,tdeg(x,y,z))[2],G);


> S:=expand(series(1/((1-x*u)*(1-y*u)*(1-z*u)),u,5));
>  L:=[op(subs(u=1, convert(S,polynom)))]; nops(L);
> with(Groebner); notdiv:={}:
> for m in L do 
>  if normalf(m,LM,plex(x,y,z))<>0 then    
>    notdiv:=notdiv union {m}: 
>  fi:
> od; print("number of solutions",nops(notdiv));
> with(Ore_algebra): A:=poly_algebra(x,y,z);
> with(Groebner): T:=termorder(A,plex(x,y,z)):
> sort(L,(t1,t2)->testorder(t1,t2,T));
> T:=termorder(A,tdeg(x,y,z));
> sort(L,(t1,t2)->testorder(t1,t2,T));
> T:=termorder(A,wdeg([3,2,1],[x,y,z])):
> sort(L,(t1,t2)->testorder(t1,t2,T));
> with(linalg): with(Groebner);

> T:=termorder(A,'matrix'([[1,1,1],[1,0,0],[0,1,0]],[x,y,z])):
> Lgrlex2:=sort(L, (t1,t2) ->testorder(t1,t2,T));