I have a system of algebraic equations. For example: $$ \left\{ \begin{aligned} x^2 y + 2 x y + 2 &= 0,\\ y^3 + 2 x + 1 &= 0. \end{aligned} \right. $$ (In my problem the system contains 16 equations of degrees 1–3).
I need to know the number of complex roots (the system in question has a finite number of roots).
I believe, the number of roots can be determined using GroebnerBasis. But what is the most elegant way to do it?
The simplest method, that will usually be correct, is to just use NSolve.
polys = {x^2 y + 2 x y + 2, y^3 + 2 x + 1}; Length[NSolve[polys]] (* Out[421]= 7 *) Below is from code I have used for solution counting and related (it is related to code inside NSolve, to answer the question of "Why not just use NSolve?"). Some is also used in this MSE response. The underlying method is described in Gröbner Bases: A Computational Approach to Commutative Algebra, by T. Becker and V. Weispfenning (with H. Kredel).
normalSet[elist_List, maxlist_, nvars_Integer] /; Length[elist] == nvars := Module[{jj, indices, iterators}, indices = Array[jj, {nvars}]; iterators = Table[{indices[[j]], 0, maxlist[[j]] - 1}, {j, nvars}]; Flatten[Table[indices, Evaluate[Apply[Sequence, iterators]]], nvars - 1]] getRoofTuples[elist_, height_] := Module[{len = Length[elist], expvec, expvecs, tt, rest}, expvecs = tt[]; Do[expvec = elist[[jj]]; If[expvec[[height]] =!= 0 && (height === 1 || Max[Take[expvec, height - 1]] === 0), expvecs = tt[expvec, expvecs]], {jj, len}]; expvecs = Apply[List, Flatten[expvecs, Infinity, tt]]; rest = Complement[elist, expvecs]; expvecs = Map[Drop[#, height - 1] &, expvecs]; {expvecs, rest}] isUnder[t1_, t2_] := Apply[And, Thread[t1 <= t2]] hitRoof[tup_, roofs_] := Module[{jj, len = Length[roofs]}, For[jj = 1, jj <= len, jj++, If[isUnder[roofs[[jj]], tup], Return[True]]]; False] normalSet[elist_List, maxes_, nvars_Integer] := Module[{tuples, ntups, tuple, ntup, mixedexpons, tt}, mixedexpons = Select[elist, (Max[Apply[Sequence, #]] != Apply[Plus, #]) &]; tuples = Table[{j}, {j, 0, maxes[[nvars]] - 1}]; Do[{roofs, mixedexpons} = getRoofTuples[mixedexpons, jj]; ntups = tt[]; Do[tuple = tuples[[kk]]; ntups = tt[Prepend[tuple, 0], ntups]; Do[ntup = Prepend[tuple, mm]; If[hitRoof[ntup, roofs], Break[], ntups = tt[ntup, ntups]], {mm, maxes[[jj]] - 1}], {kk, Length[tuples]}]; tuples = Flatten[ntups, Infinity, tt], {jj, nvars - 1, 1, -1}]; Sort[Apply[List, tuples]]] solutionSetSize[polys_, vars_] := Catch[Module[{lv = Length[vars], leadvecs, dtl, maxes, nset, pureterms}, gb = GroebnerBasis[polys, vars, MonomialOrder -> DegreeReverseLexicographic]; dtl = First[ GroebnerBasis`DistributedTermsList[gb, vars, MonomialOrder -> DegreeReverseLexicographic]]; leadvecs = Map[#[[1, 1]] &, dtl]; pureterms = Select[leadvecs, Count[#, a_ /; a =!= 0] === 1 &]; If[Length[pureterms] < lv, Throw[Infinity]]; maxes = Apply[Plus, pureterms]; nset = normalSet[leadvecs, maxes, lv]; Throw[Length[nset]] ]] Back to the example:
solutionSetSize[polys, {x, y}] (* Out[423]= 7 *) 
normalSet[elist_List, maxes_, nvars_Integer], what does normalSet do, what is a typical elist, and what sort of thing is maxes? Can you give an example? TIA $\endgroup$ solutionSetSize counts solutions by computing the normal set (with respect to a term order). I suspect that's what was wanted for the "standard monomials" in your recent question. Maybe I misinterpreted though. $\endgroup$ The documentation is very instructive.
For your particular example:
gb = GroebnerBasis[{x^2 y + 2 x y + 2, y^3 + 2 x + 1}, {x, y}] Finding roots:
sol = N@Solve[gb[[1]] == 0, y] =>{{y -> -1.29819}, {y -> -0.871821 - 1.29622 I}, {y -> -0.871821 + 1.29622 I}, {y -> 0.274176 - 1.20102 I}, {y -> 0.274176 + 1.20102 I}, {y -> 1.24674 - 0.331078 I}, {y -> 1.24674 + 0.331078 I}}
sol2 = Solve[gb[[2]] == 0, x] /. sol => {{{x -> 0.593927}}, {{x -> -2.36593 + 0.388879 I}}, {{x -> -2.36593 - 0.388879 I}}, {{x -> 0.0829231 - 0.730783 I}}, {{x -> 0.0829231 + 0.730783 I}}, {{x -> -1.26396 + 0.753779 I}}, {{x -> -1.26396 - 0.753779 I}}}
Just for fun showing real root:
ContourPlot[ Evaluate@Thread[{x^2 y + 2 x y + 2, y^3 + 2 x + 1} == 0], {x, -10, 10}, {y, -10, 10}, Epilog -> {PointSize[0.02], Red, Point[{{First[x /. First[sol2]], First[y /. sol]}}]}] 
{x^2 y + 2 x y + 2 == 0, y^3 + 2 x + 1 == 0}But the question is not about solving this particular system. $\endgroup$