# 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 '14 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 '14 at 11:43
@bcp Getting help is sometimes about making it easy for potential answerers to verify their ideas on your illustrative example. – Michael E2 Jul 2 '14 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 *)

-

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 '14 at 13:48