Computing the genus of an algebraic curve

For computing the genus of a plane algebraic curve implicitly defined by a squarefree polynomial $f(x,y)$ there are different softwares available in the literature.

Remark: I assume you are interested in computing the genus when the coefficients of the defining polynomial of the curve are either integers or rationals right? This is the case of exact data. For instance, $f(x,y)=x^2-y^3$ represents an exact data polynomial, but $f'(x,y)=1.00 \, x^2-1.01 \, y^3$, with $n=0.01$ the noise in the coefficients is an inexact data polynomial. Basically, the algorithms/their implementations for computing the genus are divided into two main classes: the algorithms for exact data and the algorithm for exact and inexact data.

Alternatives: In the case of exact data you can use some libraries, I find the following the most practical:

1. algcurves in Maple. See the author's instructions on how to compute the genus.

2. normal.lib package in Singular

In Mathematica there is no implementation (at present time from my knowledge) of an algorithm for computing the genus of plane algebraic curve. There are several reasons for this, if you want to know more about these details, just let me know.

We also have an implementation for computing the genus of a plane algebraic curve at, please see here.

Roughly, the implementation is in C++ (see the instructions for more details). Our library computes the genus of curves defined by polynomials with exact data and with inexact data. We did not write our library in Mathematica for different "strong" reasons, so if you really want to use Mathematica for computing the genus I am afraid I do not have a good news for you. But, for the exact case, the algcurves package is very good, and for the exact and inexact case (or simply the exact case) you can use our library. I hope this is not too confusing for you.

If you need more help/information, please just let me know.

Note: Here are some reasons for which we did not implement at first our algorithm for genus computation in Mathematica (Artes, thanks for the question on this issue, I was kind of hoping to discuss a bit more about this).

For computing the genus, we need an algorithm for solving the following decisive problem: for a real plane algebraic curve C defined by the polynomial f(x,y)=0, we need to compute a graph G=(V,E), where V is a set of points in the 2-dimensional Euclidean plane together with their Euclidean coordinates and E is a set of edges connecting them. In addition, the graph G has the special property that it can be continuously deformed into the curve C. Intuitively, the computed graph G is a piecewise linear approximation of the differential curve C. At present (from my knowledge) in Mathematica there are no algorithms for computing the graph G for an input curve C. Still, at present we are working (I am actually very interested in this problem) on developing a new algorithm for computing such a special graph G for a real plane (and space) real algebraic curve. And I am working on the implementation of the algorithm in Mathematica. After this algorithm will be available in Mathematica, we can have an implementation for the genus of a plane algebraic curve in Mathematica as well.

I will show a method that is conjectural, though i believe it is correct. It differs from the more common approach of using (quadratic) birational transformations to force singularities to be double points. More a detailed approach to that, see Madelina's response. Also I cover the exact case, although at least some of this could be adopted to the case of approximate coefficients.

Recall the genus formula $ g=\binom{d-1}{2}-\sum _ {m_p\in S} \binom{m_p}{2} $ where $S$ is the set of singular points on the curve, and $m_p$ is the multiplicity of point $p$. There is a catch of sorts: the multiplicity is not in general the same as what one obtains from solving the appropriate polynomial system to find the singularities. That only holds for what are called "ordinary" multiple points. So the game plan is to find the singularities and the correct multiplicities. Here is a method that I believe works although I do not have a proof for this. The idea is as follows.

(1) Start with a polynomial $p(x,y)$ that implicitly gives a planar curve.

(2) Perform a random linear change of coordinates to move the center and also to change the directions of $x$ and $y$. The change of center is for numeric purposes: it is sometimes difficult to assess multiplicity of numeric solutions at the origin. The change of directions is to be sure that no coordinate has "accidental" multiplicity, that is, apepars as a solution value in two or more distinct solutions of a certain system of equations to be discussed.

(3) Homogenize the new polynomial in a new variable $z$.

(4) Project the homogenized polynomial onto a random affine hyperplane. The reason for this is that "badness' might be happening in projective 2-space "at infinity", that is, off the $z=1$ plane where the polynomial shows its $(x,y)$ persona. By projecting onto a random affine hyperplane we make sure any such badness will appear in our new affine coordinates.

For a simple example of why this is at all important, consider $p(x,y) = y-x^4$. It is obviously parametrizable (by $(x,x^4)$), hence has genus 0. Also it has no (finite) singularities because there is no place where $D_y p(x,y)$ vanishes. The genus formula would therefore claim a genus of $(d-1)(d-2)/2$ where $d=degree(p)$, hence 3.

(5) Set up the polynomial system that locates singular points: ${p(x,y), D_x p (x, y), D_y p (x, y)}$ Now compute a Groebner basis for this system. While we're at it, factor the polynomials in this basis. It might help to see what is happening, at least for the standard uglyissimo example to follow.

(6) Not needed, but we'll use NSolve to find the solutions to the singularity equations. This will also show multiplicity of those solutions. Which, alas, is not quite what we need in some cases.

(7) We've randomized enough to satisfy any realistic needs, so we now set up a matrix whose characteristic polynomial gives exactly the values of $y$ at the singularities (ordinarily I would take a random linear combination of the variables, but prior randomization has obviated the need for that). This matrix is referred to as a generalized characteristic matrix, if I recall correctly. Assuming we computed a basis with respect to lexicographic term order, with variables ordered so that $x>y$, the characteristic polynomial of the matrix will be a close cousin to the first polynomial in our basis (roots will agree, degrees may not).

I will remark that the code to do this comes from implementation of NSolve. It is a bit long but not terribly relevant to the main ideas here. There is an existing body of literature on solving polynomial systems by eigensystems, and indeed it generalizes the method of finding polynomial roots numerically by computing eigenvalues of the matrix formed from the characteristic polynomial.

(8) Compute the JordanDecomposition of this matrix. Use that to deduce the factored minimal polynomial of the variable in question. The degrees of each root of this minimal polynomial will give the correct multiplicities to use in the genus formula. This is the part for which I do not have a proof, but I am fairly sure it works.

Computationally, this last might be the most troublesome step since it involves linear algebra possibly in extension fields of the rationals. There may be ways to avoid this using only the rational canonical form of the generalized characteristic matrix.

With all that said, here is an example. It is known to have genus zero. As I found it in homogenized form, I'll show it that way even though we will straight off dehomogenize, do our transformations, and only then rehomogenize.

In[72]:= vars = {x, y};
allvars = Union[Append[vars, z]];
SeedRandom[123456];
center = RandomInteger[{1, 20}, Length[vars]]
det = 0;
While[det == 0, 
  det = Det[transformmat = RandomInteger[{-5, 5}, {2, 2}]]];
transformmat
det

Out[75]= {16, 3}

Out[78]= {{2, 1}, {-4, -5}}

Out[79]= -6

So we have a viable transformation to work with below.

poly = x^5 + 10*x^4*y + 20*x^3*y^2 + 130*x^2*y^3 - 20*x*y^4 + 
   20*y^5 - 2*x^4*z - 40*x^3*y*z - 150*x^2*y^2*z - 90*x*y^3*z - 
   40*y^4*z + x^3*z^2 + 30*x^2*y*z^2 + 110*x*y^2*z^2 + 20*y^3*z^2;
poly2 = Expand[(poly /. z -> 1) /. 
    Thread[vars -> (transformmat.vars - center)]];
cp = ContourPlot[poly /. z -> 1, {x, -3, 3}, {y, -1, 2.5}, 
  Contours -> {0}, ContourShading -> False, Frame -> False, 
  PlotPoints -> 80, ImageSize -> 300]

We will show a picture of the zero contours of the two first derivatives as well, zooming in on the origin. One can see the funky intersection behavior at the origin.

Here I perturb the two derivative curves and zoom in further. We see a few salient features. One (perhaps already observed) is that the polynomial itself seems to have three mutually nontangential crossings. I believe this indicates a triple point singularity, although my curvology is, frankly, quite weak. Another is that the intersections with derivatives appears to have some tangential crossings. One can corroborate this by seeing that the second derivative zero contours also intersect there. I do not know what is the significance of this. A third observation is that the multiplicity of the intersections with the three curves might be four, since we have sets of pairwise intersections that appear to merge, and I'm guessing pairs come together as one when the derivatives are left unperturbed. This is all of course just speculation, based on graphical information.

We proceed to do the randomization, homogenization, and a few other izations.

In[1063]:= hpoly = homogenize[poly2, vars, z];
SeedRandom[111222];
randomplane = (RandomInteger[{10, 100}, 3]*
      RandomChoice[{-1, 1}, 3]).allvars - 100;
newpoly = First[GroebnerBasis[{hpoly, randomplane}, vars, z]];
polys = {newpoly, D[newpoly, x], D[newpoly, y]};
gb = GroebnerBasis[polys, vars];
Factor[gb]

Out[1069]= {(-925 + 1182 y) (-1200 + 1519 y) (-7000 + 
    8907 y)^3 (-7600 + 9663 y), (-7000 + 
    8907 y) (892380962668062655597889053350846133200000000 + 
    1114037523664535012054078581248000000 x - 
    4540867588642210093905227481286786853197000000 y + 
    8664769369921590249333220599218979749558700000 y^2 - 
    7348370038737755798344798868620889302139772225 y^3 + 
    2336978157116247652016918398962438344469882294 y^4), \
-283829738542292200993193596906934417400491846960000000000 + 
  46249953002279392958759819695776759603200000000 x + 
  24816164541644732113474320122306240830464000000 x^2 + 
  1805426407890987718981716357721013439583910821161000000000 y - 
  4593679607761292838195012880861431289086833972537643000000 y^2 + 
  5844005854319148261676988351692967316344655988400545875000 y^3 - 
  3717317126255736588220159220821056995508100012116655092675 y^4 + 
  945818051707434857990222520578342245206067344137050058178 y^5}

One will notice that one of the factors in the first term appears to degree higher than one, and also appears as a factor in the next polynomial. This is, in the technical parlance, "bad news". It means the characteristic matrix was derogatory; we will see that explicitly in the Jordan form below. First we show the numeric values of the singularities, along with their algebraic multiplicity.

In[1070]:= Gather[NSolve[gb == 0, WorkingPrecision -> 400]] // N

Out[1070]= {{{x -> -0.930626, y -> 0.782572}}, {{x -> -0.921659, 
   y -> 0.789993}}, {{x -> -0.921039, 
   y -> 0.786505}}, {{x -> -0.931851, y -> 0.785899}, {x -> -0.931851,
    y -> 0.785899}, {x -> -0.931851, y -> 0.785899}, {x -> -0.931851, 
   y -> 0.785899}}}

Here I form the generalized characteristic matrix. The somewhat lengthy code for all this, and for the homogenization used above, is at the end of this note.

ordr = Lexicographic;
len = Length[vars];
{heads, tails, dim} = breakApart[gb, vars, ordr];
elist = Map[First, heads];
maxexpons = Map[Max, Transpose[elist]];
headmonoms = headsTable[elist, maxexpons, len];
hlen = Length[headmonoms];

(* Had we not earlier randomized everything, here we would use
varvec = getRandomVector[len]; *)
varvec = {0, 1};
endomat1 = 
  getEndoMatrix[headmonoms, varvec, hlen, vars, gb, Infinity, ordr];

The characteristic polynomial:

In[1199]:= Factor[CharacteristicPolynomial[endomat1, y]]

Out[1199]= -(((-925 + 1182 y) (-1200 + 1519 y) (-7000 + 
    8907 y)^4 (-7600 + 9663 y))/109197586392254572178903454)

So the "bad" root has degree 4 here. Let's see what the Jordan normal form of the characteristic matrix has to say.

In[1085]:= JordanDecomposition[endomat1][[2]]

Out[1085]= {{-(2274447700/1145360037), 0, 0, 0, 0, 0, 
  0}, {0, -(2274447700/1145360037), 1, 0, 0, 0, 0}, {0, 
  0, -(2274447700/1145360037), 1, 0, 0, 0}, {0, 0, 
  0, -(2274447700/1145360037), 0, 0, 0}, {0, 0, 0, 
  0, -(301084525/151994562), 0, 0}, {0, 0, 0, 0, 
  0, -(55159800/27904247), 0}, {0, 0, 0, 0, 0, 
  0, -(2450940100/1242574833)}}

One will observe that the $(-7000 + 8907 y)^4$ factor in the Characteristic Polynomial will only be of degree 3 in the Minimal Polynomial. This follows from the fact that the corresponding eigenvalue splits into two Jordan blocks of length 3 and 1 respectively.

From this follows the contention that the corresponding root is a triple point singularity, hence we should subtract $\binom{3}{2}=3$ for this singular point. As the remaining three singularities are of multiplicity 1, we subtract 1 for each of them. Since the degree $d=5$ the first term is $\binom{4}{2}=6$. Hence the genus is $6-3-1-1-1=0$.

This approach raises several questions, the most obvious being, is it really correct? Another is whether the Groebner basis might be used in lieu of the eventual Jordan form, to deduce the "correct" multiplicity to use in the genus formula. Yet another is to understand, from a geometric viewpoint, the cause of the derogatoriness of the matrix in question. I do not have answers so I'll leave off here.

As promised, here is the full code used above.

DTL = GroebnerBasis`DistributedTermsList;

breakApart[basis_, pvars_, order_] := 
 Module[{exponvecs, heads, rest, dim, polys},
  polys = 
   First[DTL[basis, pvars, Sort -> False, MonomialOrder -> order, 
     CoefficientDomain -> RationalFunctions]];
  Quiet[heads = 
    Map[First, polys] /. HoldPattern[First[{}]] :> Sequence[]];
  Quiet[rest = 
    Map[Rest, polys] /. HoldPattern[Rest[{}]] :> Sequence[]];
  exponvecs = Map[First, heads];
  exponvecs = exponvecs /. n_?Positive -> 1;
  exponvecs = Select[exponvecs, Apply[Plus, #] === 1 &];
  dim = ! 
    SameQ[Reverse[Sort[exponvecs]], IdentityMatrix[Length[pvars]]];
  {heads, rest, dim}]

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

headsTable[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]]]

coeffExpon[{hmonom_, coeff_}, hmonoms_] := 
 Module[{len = Length[hmonoms], pos}, 
  pos = Position[hmonoms, hmonom, {1}, 1];
  If[pos === {}, Return[$Failed]];
  {pos[[1, 1]], coeff}]

coefficientRow[monoms_, hmonoms_] := 
 Module[{clist, len = Length[hmonoms], ncoeffs, kk = 1, coeffs, 
   expons}, clist = Sort[Map[coeffExpon[#, hmonoms] &, monoms]];
  If[! FreeQ[clist, $Failed], Return[$Failed]];
  {expons, coeffs} = Transpose[clist];
  ncoeffs = Length[expons];
  Table[If[kk <= ncoeffs && j == expons[[kk]], kk++; coeffs[[kk - 1]],
     0], {j, len}]]

getEndoMatrix[headmonoms_, varvec_, hlen_, pvars_, gb_, wprec_, 
  order_] := 
 Module[{newterms, nvars = Length[pvars], headterms, crow, res}, 
  headterms = Map[pvars^# &, headmonoms];
  headterms = Map[Apply[Times, #] &, headterms];
  newterms = Sum[varvec[[k]]*pvars[[k]]*headterms, {k, nvars}];
  Quiet[newterms = 
    PolynomialReduce[newterms, gb, pvars, 
     CoefficientDomain -> 
      If[wprec =!= Infinity, InexactNumbers, RationalFunctions], 
     MonomialOrder -> order]];
  If[Head[newterms] === PolynomialReduce, Return[$Failed]];
  newterms = Map[Last, newterms];
  res = Table[newterm = newterms[[j]];
    If[PossibleZeroQ[newterm], Table[0, {hlen}],(*else*)

     newterm = 
      First[DTL[newterm, pvars, Sort -> False, 
        MonomialOrder -> order]];
     crow = coefficientRow[newterm, headmonoms];
     crow], {j, hlen}];
  res]

homogenize[poly_, vars_, new_] := 
 Module[{degfunc, totdeg, j}, 
  degfunc = Apply[Plus, Map[Function[x, Exponent[#, x] &], vars]];
  degfunc = Distribute[degfunc, Function, Plus];
  totdeg = Max[Map[degfunc, Apply[List, poly]]];
  Apply[Plus, 
   Table[poly[[j]]*new^(totdeg - degfunc[poly[[j]]]), {j, 
     Length[poly]}]]]

getRandomVector[len_] := 
 Module[{vec, den}, 
  vec = (2*RandomInteger[{0, 1}, len] - 1)*
    RandomInteger[{100000, 200000}, len];
  den = RandomInteger[{100000, 200000}];
  vec/den]

Mădălina Hodorog has software to do this.

She has some Mathematica packages too, so she might really have the last word on this.

She publishes her email address, so you can contact her directly at the Technical University Berlin Institute of Mathematics:

[email protected]

And a list of her publications.