# Is there a function in Mathematica that computes the number of solutions to polynomial systems?

Given a system of polynomial equations in $$\mathbb{C}$$-coefficients, is there a tool in Mathematica that computes the number of solutions to this system, counted with multiplicity? (We may assume that the algebraic variety defined by the system is zero-dimensional.)

Equivalently, assuming we are given an ideal $$I\subset \mathbb{C}[x_1,\dots,x_n]$$, and it is known that $$\mathbb{C}[x_1,\dots,x_n]/I$$ is a finite dimensional linear space over $$\mathbb{C}$$. Is there is a tool in Mathematica that computes the complex dimension of the quotient space?

This can be done by Maple using the NumberOfSolutions function in the PolynomialIdeals package, but I am trying to learn how to do it in Mathematica. It might be somewhere in the GroebnerBasis package, but I can't seem to find it in the online manual.

Thanks a lot!

Here is a concrete example:

If we are given two polynomials x^2-y, y^3+y+1, then the answer should be 6, because the system of equations $$x^2-y=0, y^3+y+1=0$$ has 6 solutions when counted with multiplicity.

I believe the algorithm needs to somehow compute the Groebner basis of the given ideal, but I am not familiar with the details of the algorithms and couldn't figure out how to use the Groebner package for it.

Thanks again!

• A concrete example always helps. Mar 5, 2023 at 22:14
• NSolve can be used with a certain Method setting. Mar 5, 2023 at 22:48
• Have you tried CountRoots? If you have why doesn't it help? Perhaps also RootIntervals could be helpful. Mar 6, 2023 at 0:13
• Did you see GroebnerBasis[{x^2 - y, y^3 + y + 1}, {y}] or Eliminate[{x^2 - y == 0, y^3 + y + 1 == 0}, {y}]?
– Moo
Mar 6, 2023 at 11:17
• Please edit the post to plut the example system in copy-pastable format (Mathematica InputForm works well for this). It is the expectation for a minimal working example in this forum, and saves people from retyping or having to check comments or responses to find such code. Mar 6, 2023 at 15:38

Here is code I cribbed from the NSolve method based on computing an eigensystem from a Groebner basis.

borderMonomials[basis_, pvars_, order_] := Module[
polys =
First[GroebnerBasisDistributedTermsList[basis, pvars,
Sort -> False, MonomialOrder -> order,
CoefficientDomain -> RationalFunctions]];
HoldPattern[First[{}]] :> Sequence[]];

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_] :=
Catch[Module[{jj, len = Length[roofs]},
For[jj = 1, jj <= len, jj++,
If[isUnder[roofs[[jj]], tup], Throw[True, "hR"]]];
False], "hR"]

Module[{tuples, ntups, tuple, ntup, mixedexpons, tt, roofs},
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]]]

CountSolutions[polys_, vars_] := Module[
gb = GroebnerBasis[polys, vars,
CoefficientDomain -> RationalFunctions,
MonomialOrder -> DegreeReverseLexicographic];
elist = borderMonomials[gb, vars, DegreeReverseLexicographic];
maxexpons = Map[Max, Transpose[elist]];
]


Here is the example from the post.

polys = {x^2 - y, y^3 + y + 1};
vars = Variables[polys];
CountSolutions[polys, vars]

(* Out[142]= 6 *)


Here is a more complicated example.

polys = {x^3 + 5*x*z - 2*x*y - y^2 + z^3 - 3*y^2*z,
y^3 + x^2*y + -7*x*y*2 + 3*x^2 - 5*x*y - 1,
x*z^2 + 4*x*y*z + y^2*z - 2*z + x^2*y*z - 11};
vars = Variables[polys];
CountSolutions[polys, vars]

(* Out[155]= 33 *)


When I get a chance I'll prepare and submit this to the Wolfram Function Repository.

• fyi, using the Maple command the OP linked above in the comment, it gives different answer for one same input test I did. Here is an example. polys=={x^2-y+x^9;y^8+y+1+x^4+2*y*x}; vars=Variables[polys]; CountSolutions[polys,vars] gives 6 but Maple gives $\infty$ for same polynomials. Here is screen shot !Mathematica graphics Mar 6, 2023 at 20:13
• Thaks for that example @Nasser. Yes, infinity is correct. I should have noted that my code explicitly expects/works on systems with finitely many solutions (this was stipulated in the original post). It is not hard to adjust to be more general and indeed NSolve code does have the necessary checks in place. Mar 6, 2023 at 20:47
• Thank you so much for the code! It seems that Maple and your code give different answers for the following system: polys = {a^5 + 30*a^3*b + 40*a^2*c + 89*ab^2 + 88*bc, a^6 + 55*a^4*b + 80*a^3*c + 439*a^2*b^2 + 688*abc + 225*b^3 + 160*c^2, a^7 + 91*a^5*b + 140*a^4*c + 1519*a^3*b^2 + 2968*a^2*bc + 3429*ab^3 + 1120*a*c^2 + 3708*b^2*c} Your code gives 146 but Maple gives 35. Sorry for this complicated example, it is actually produced by some other computations and for this example, we do know that the answer is 35. Maybe there is just some minor bug. Thanks again! Mar 7, 2023 at 3:37
• @BoyuZhang I get 35 from Daniel's code if I put a space in between the letters in ab, abc, etc. Mar 7, 2023 at 15:08
• Daniel, I don't know if you want your WFR function to handle the following, but it gives and error and then the correct answer: CountSolutions[{x^3 + Sqrt[y^4] - 1, x^2 + y^2 - 1}, {x, y}]` Mar 7, 2023 at 15:09