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

Question

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!

Answer 1

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

borderMonomials[basis_, pvars_, order_] := Module[
  {polys, heads},
  polys = 
   First[GroebnerBasis`DistributedTermsList[basis, pvars, 
     Sort -> False, MonomialOrder -> order, 
     CoefficientDomain -> RationalFunctions]];
  Quiet[heads = Map[First, polys] /.
     HoldPattern[First[{}]] :> Sequence[]];
  Map[First, heads]]

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

headsTable[elist_List, maxes_, nvars_Integer] := 
 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, elist, maxexpons, headmonoms},
  gb = GroebnerBasis[polys, vars, 
    CoefficientDomain -> RationalFunctions, 
    MonomialOrder -> DegreeReverseLexicographic];
  elist = borderMonomials[gb, vars, DegreeReverseLexicographic];
  maxexpons = Map[Max, Transpose[elist]];
  headmonoms = headsTable[elist, maxexpons, Length[vars]];
  Length[headmonoms]
  ]

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.