Sextic Classifier

Following is a software in Mathematica that determines the topological type of a given plane sextic. The input is a homogeneous polynomial f of degree six in three variables x,y,z and the output is the topological type. Note that if your input is a singular curve output would be “Singular”, otherwise it is of the forms: emptyset, {n}, {n,1}, {{n,1},m}, and {1,1,1} where m and n are integers.

Here {n} means n ovals with no nesting, {n,1} is n ovals sitting outside eachother and all nested in another oval, {{n, 1}, m} means there is an oval with n ovals inside it and outside each other and m ovals outside it and outside each other and finally {1,1,1} is the hyperbolic type where we have an oval nested in another oval and together, they are nested in the third oval.

CheckEqual[g_] :=
Not[Length[g] > 2] && (g[[0]][1, 1] && Not[g[[0]][0, 1]] &&
Not[g[[0]][1, 0]])

Fiber[cyl_, i_] :=
Module[{res},
If[Length[cyl[[i, 2, 1]]] == 0, res = 1, res = Length[cyl[[i, 2]]]];
res]

Base[cyl_] := Table[i, {i, 1, Length[cyl]}]

Vertices[cyl_] :=
Module[{ram, bra, i, j}, ram = {}; bra = Base[cyl];
For[i = 1, i <= Length[bra], i++,
For[j = 1, j <= Fiber[cyl, bra[[i]]], j++,
ram = Append[ram, {bra[[i]], j}]]]; ram]

Branch[cyl_] :=
Module[{bra, i}, bra = {};
For[i = 1, i <= Length[cyl], i++,
If[CheckEqual[cyl[[i, 1]]], bra = Append[bra, i]]]; bra]

RamNumber[cyl_, i_] :=
Module[{che, j},
If[Fiber[cyl, i] == 1,
che = 1, {che = Fiber[cyl, i];
For[j = 1, j < Fiber[cyl, i], j++,
If[cyl[[i, 2, j, 2, 2]] < cyl[[i, 2, j + 1, 2, 2]] - 1,
che = j]]};]; che]

EdgeNBtoBMore[cyl_, i_] :=
Module[{ram, j, G, met}, met = 0; G = {};
ram = RamNumber[cyl, i + 1];
For[j = 1, j <= Fiber[cyl, i], j++, If[j == ram, met = 1];
G = Append[G, {i, j} <-> {i + 1, j + met}]]; G]

EdgeNBtoBLess[cyl_, i_] :=
Module[{ram, j, G, met}, met = 0; G = {};
ram = RamNumber[cyl, i + 1];
For[j = 1, j <= Fiber[cyl, i], j++,
G = Append[G, {i, j} <-> {i + 1, j - met}];
If[j == ram, met = 1];]; G]

EdgeBtoNBMore[cyl_, i_] :=
Module[{ram, j, G, met}, met = 0; G = {}; ram = RamNumber[cyl, i];
For[j = 1, j <= Fiber[cyl, i + 1], j++,
G = Append[G, {i, j - met} <-> {i + 1, j}]; If[j == ram, met = 1]];
G]

EdgeBtoNBLess[cyl_, i_] :=
Module[{ram, j, G, met}, met = 0; G = {}; ram = RamNumber[cyl, i];
For[j = 1, j <= Fiber[cyl, i + 1], j++, If[j == ram, met = 1];
G = Append[G, {i, j + met} <-> {i + 1, j}];]; G]

Edges[cyl_] :=
Module[{i, t, G, j}, G = {}; For[i = 1, i < Length[cyl], i++,
If[CheckEqual[cyl[[i, 1]]],
If[CheckEqual[cyl[[i + 1, 1]]], t = t,
If[Fiber[cyl, i] < Fiber[cyl, i + 1],
G = Union[G, EdgeBtoNBMore[cyl, i]],
G = Union[G, EdgeBtoNBLess[cyl, i]]]],
If[Fiber[cyl, i] < Fiber[cyl, i + 1],
G = Union[G, EdgeNBtoBMore[cyl, i]],
G = Union[G, EdgeNBtoBLess[cyl, i]]]]]; i = Length[cyl];
If[CheckEqual[cyl[[i, 1]]], t = t,
For[j = 1, j <= Fiber[cyl, i], j++,
G = Append[G, {i, j} <-> {1, Fiber[cyl, i] + 1 - j}]];]; G]

EdgesAffine[cyl_] :=
Module[{i, t, G, j}, G = {}; For[i = 1, i < Length[cyl], i++,
If[CheckEqual[cyl[[i, 1]]],
If[CheckEqual[cyl[[i + 1, 1]]], t = t,
If[Fiber[cyl, i] < Fiber[cyl, i + 1],
G = Union[G, EdgeBtoNBMore[cyl, i]],
G = Union[G, EdgeBtoNBLess[cyl, i]]]],
If[Fiber[cyl, i] < Fiber[cyl, i + 1],
G = Union[G, EdgeNBtoBMore[cyl, i]],
G = Union[G, EdgeNBtoBLess[cyl, i]]]]]; i = Length[cyl]; G]

Comps[cyl_] := Length[ConnectedComponents[Graph[Edges[cyl]]]]

FiberOval[A_, k_] :=
Module[{i, count}, count = 0;
For[i = 1, i <= Length[A], i++,
If[A[[i, 1]] == k, count = count + 1];]; count];

OvalTypeAffine[A_, cyl_] :=
Module[{res, i}, res = 1;
For[i = 1, i <= Length[cyl], i++, If[FiberOval[A, i] < 1, res = 0]];
res]

(*OvalType[A_,cyl_]:=Module[{B,j,V,type1,type2},res=1;type1=0;type2=0;\
B=ConnectedComponents[Graph[EdgesAffine[cyl]]];For[j=1,j\[LessEqual]\
Length[B],j++,If[OvalTypeAffine[B[[j]],cyl]\[Equal]1,V=B[[j,1]];If[\
MemberQ[A,V],type1=type1+1],V=B[[j,1]];If[MemberQ[A,V],type2=type2+1]]\
];If[CheckEqual[cyl[[1,1]]],res=0];
If[type2>0,res=0];
res]*)

RandVertices[A_, len_] :=
Module[{list, i}, list = {};
For[i = 1, i <= Length[A], i++,
If[A[[i, 1]] == 1 || A[[i, 1]] == len,
list = Append[list, A[[i]]]]]; list];

RandVerticesOther[A_, len_, rand1_] :=
Module[{list, i, res}, list = {};
For[i = 1, i <= Length[A], i++,
If[A[[i, 1]] == 1 || A[[i, 1]] == len,
list = Append[list, A[[i]]]]];
If[list[[1]] == rand1, res = list[[2]], res = list[[1]]]; res];

FindComp[B_, vert_] :=
Module[{i, res},
For[i = 1, i <= Length[B], i++, If[MemberQ[B[[i]], vert], res = i]];
res]

OvalType[A_, cyl_] :=
Module[{B, j, V, len, findstart, angle, direction, finish, start,
startvert, li},
B = ConnectedComponents[Graph[EdgesAffine[cyl]]];
findstart = 0; j = 1;
While[findstart == 0,
V = B[[j, 1]];
If[MemberQ[A, V], findstart = 1; start = j;, j = j + 1];
];
If[(Length[A]) == (Length[B[[j]]]), res = 0,
{direction = 1;
angle = 0;
finish = 0;
len = Length[cyl];
j = start;
li = RandVertices[B[[j]], len];
startvert = li[[1]];
If[OvalTypeAffine[B[[start]], cyl] == 1 && startvert[[1]] > 1,
startvert = li[[2]]];
While[finish == 0,
If[OvalTypeAffine[B[[j]], cyl] == 1, angle = angle + direction,
direction = -direction];
startvert = RandVerticesOther[B[[j]], len, startvert];
startvert = {1 + len - startvert[[1]],
Fiber[cyl, len] + 1 - startvert[[2]]};
j = FindComp[B, startvert];
If[j == start, finish = 1];];
If[angle == 2, res = 1, res = 0];
If[CheckEqual[cyl[[1, 1]]], res = 0]}];
res]

Nested[A_, B_, cyl_] :=
Module[{count, i, bra, V, j}, bra = Branch[cyl];
For[i = 1, i <= Length[A], i++,
If[Not[MemberQ[bra, A[[i, 1]]]], j = i]];
V = A[[j]]; count = 0;
For[i = 1, i <= Length[B], i++,
If[V[[1]] == B[[i, 1]] && V[[2]] < B[[i, 2]], count = count + 1]];
OddQ[count + OvalType[B, cyl]]]

TopType[f_] :=
Module[{edges, G, comps, i, j, type, count, sol, f0, cyl},
cyl = CylindricalDecomposition[{(f /. {z -> 1}) == 0}, {x, y}];
If[NoBranch[f], f0 = f /. {x -> 0, z -> 1};
sol = Solve[f0 == 0, Reals]; comps = Length[sol]/2;
If[comps == 0, type = "Empty"]; If[comps == 1, type = {1}];
If[comps == 2, type = {{1, 1}, 0}];
If[comps == 3, type = {1, 1, 1}]; If[Sing[f], type = "Singular"],
If[Length[cyl] == 0, type = "Empty", edges = Edges[cyl];
G = Graph[edges]; comps = ConnectedComponents[G]; count = 0;
For[i = 1, i <= Length[comps], i++,
For[j = 1, j <= Length[comps], j++,
If[i != j && Nested[comps[[i]], comps[[j]], cyl],
count = count + 1]]]; If[count == 0, type = {Length[comps]}];
If[count > 0, type = {{count, 1}, Length[comps] - count - 1}];
If[Length[comps] == 3 && count == 3, type = {1, 1, 1}];
If[BadCoords[f] || TooNice[cyl],
type = "Please change the coordinates"];
If[Sing[f], type = "Singular"];];]; type]

Sing[f_] := (Length[
Union[Solve[{D[f, x], D[f, y], D[f, z], x} == {0, 0, 0, 1}],
Solve[{D[f, x], D[f, y], D[f, z], y} == {0, 0, 0, 1}],
Solve[{D[f, x], D[f, y], D[f, z], z} == {0, 0, 0, 1}]]] > 0)

Unbounded[
f_] := (Length[
Union[Solve[{f, x, z} == {0, 1, 0}, Reals],
Solve[{f, y, z} == {0, 1, 0}, Reals]]]) > 0

CenterofProjection[f_] := ((f /. {x -> 0, z -> 0, y -> 1}) == 0)

DoubleBranch[f_] :=
Module[{disc, fa}, fa = f /. {z -> 1};
disc = GroebnerBasis[{fa, D[fa, y]}, {x}, {y}][[1]];
Exponent[disc, x] < 30]

NoBranch[f_] :=
Module[{fa, sol, k}, fa = f /. {z -> 1};
sol = Solve[{fa, D[fa, y]} == 0, Reals]; Length[sol] == 0]

TooNice[cyl_] :=
Module[{i, j, res, t}, res = 0;
For[i = 1, i <= Length[cyl], i++,
If[Fiber[cyl, i] == 1, t = 0,
For[j = 1, j <= Length[cyl[[i, 2]]], j++,
If[Depth[cyl[[i, 2, j, 2]]] <= 1, res = 1]]]]; res == 1]

BadCoords[f_] := DoubleBranch[f] || CenterofProjection[f]

MakeGoodCoords[f_] :=
Module[{g, A, h}, g = f;
If[Not[Sing[f]],
While[BadCoords[g], A = RandomInteger[{-10, 10}, {3, 3}];
If[Det[A] != 0, h = A.{x, y, z};
g = g /. {x -> h[[1]], y -> h[[2]], z -> h[[3]]}]];]; g]

f = 19 x^6 - 20 x^4 y^2 - 20 x^2 y^4 + 19 y^6 - 20 x^4 z^2 +
60 x^2 y^2 z^2 - 20 y^4 z^2 - 20 x^2 z^4 - 20 y^2 z^4 + 19 z^6;

f = f /. {x -> x + 2*y};

TopType[f]

{10}