# Number of complex roots for the system of polynomials

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?

• Your question would get more attention if you provide Mathematica code instead of $\LaTeX$ :)
– Öskå
Jul 1, 2014 at 10:50
• {x^2 y + 2 x y + 2 == 0, y^3 + 2 x + 1 == 0} But the question is not about solving this particular system.
– bcp
Jul 1, 2014 at 11:43
• @bcp Getting help is sometimes about making it easy for potential answerers to verify their ideas on your illustrative example. Jul 2, 2014 at 0:49

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[
GroebnerBasisDistributedTermsList[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];
Throw[Length[nset]]
]]


Back to the example:

solutionSetSize[polys, {x, y}]

(* Out[423]= 7 *)

• In 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
– atat
Apr 23, 2023 at 1:37
• @atat Given a list of polynomials and variables, 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. Apr 23, 2023 at 1:44

The documentation is very instructive.

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]}}]}]
`

• So the question is whether I can determine the number of roots without actually finding all of them. (e.g. if I had only one single-variable polynomial equation, I would take its degree as the answer)
– bcp
Jul 1, 2014 at 13:48