(* Mathematica computation of the nine inflection points of a general ternary cubic in radicals. *)

Clear[x,y,z,c300,c210,c120,c030,c201,c111,c021,c102,c012,c003];

f=c300 x^3+3c210 x^2 y + 3c201 x^2 z + 3c120 x y^2 + 6c111 x y z + 3c102 x z^2 + c030 y^3 + 3c021 y^2 z + 3c012 y z^2 + c003 z^3;

hessian = Det[D[f,{{x,y,z},2}]]/216;

(* We have left the following formulas for the degree 4 and 6 invariants S and T unexpanded so as to save space in the resulting formulae for the inflection points. 

S = -1*(c300 c120 c021 c003 -  c300 c120 c012^2 -  c300 c111 c030 c003 +  c300 c111 c021 c012 +  c300 c102 c030 c012 -  c300 c102 c021^2 -  c210^2 c021 c003+  c210^2 c012^2 +  c210 c201 c030 c003 -  c210 c201 c021 c012 +  c210 c120 c111 c003-  c210 c120 c102 c012 - 2 c210 c111^2 c012 + 3 c210 c111 c102 c021 -  
c210 c102^2 c030-  c201^2 c030 c012 +  c201^2 c021^2 -  c201 c120^2 c003 + 3 c201 c120 c111 c012-  c201 c120 c102 c021 - 2 c201 c111^2 c021 + 
 c201 c111 c102 c030 +  c120^2 c102^2- 2 c120 c111^2 c102 +  c111^4);
T=c300^2 c030^2 c003^2-6 c300^2 c030 c021 c012 c003+4 c300^2 c030 c012^3+4 c300^2 c021^3 c003-3 c300^2 c021^2 c012^2-6 c300 c210 c120 c030 c003^2+18 c300 c210 c120 c021 c012 c003-12 c300 c210 c120 c012^3+12 c300 c210 c111 c030 c012 c003-24 c300 c210 c111 c021^2 c003+12 c300 c210 c111 c021 c012^2+6 c300 c210 c102 c030 c021 c003-12 c300 c210 c102 c030 c012^2+6 c300 c210 c102 c021^2 c012+6 c300 c201 c120 c030 c012 c003-12 c300 c201 c120 c021^2 c003+6 c300 c201 c120 c021 c012^2+12 c300 c201 c111 c030 c021 c003-24 c300 c201 c111 c030 c012^2+12 c300 c201 c111 c021^2 c012-6 c300 c201 c102 c030^2 c003+18 c300 c201 c102 c030 c021 c012-12 c300 c201 c102 c021^3+4 c300 c120^3 c003^2-24 c300 c120^2 c111 c012 c003-12 c300 c120^2 c102 c021 c003+24 c300 c120^2 c102 c012^2+36 c300 c120 c111^2 c021 c003+12 c300 c120 c111^2 c012^2+12 c300 c120 c111 c102 c030 c003-60 c300 c120 c111 c102 c021 c012-12 c300 c120 c102^2 c030 c012+24 c300 c120 c102^2 c021^2-20 c300 c111^3 c030 c003-12 c300 c111^3 c021 c012+36 c300 c111^2 c102 c030 c012+12 c300 c111^2 c102 c021^2-24 c300 c111 c102^2 c030 c021+4 c300 c102^3 c030^2+4 c210^3 c030 c003^2-12 c210^3 c021 c012 c003+8 c210^3 c012^3-12 c210^2 c201 c030 c012 c003+24 c210^2 c201 c021^2 c003-12 c210^2 c201 c021 c012^2-3 c210^2 c120^2 c003^2+12 c210^2 c120 c111 c012 c003-24 c120^2 c111^2 c102^2+24 c120 c111^4 c102+6 c210^2 c120 c102 c021 c003-12 c210^2 c120 c102 c012^2+12 c210^2 c111^2 c021 c003-24 c210^2 c111^2 c012^2-24 c210^2 c111 c102 c030 c003-27 c210^2 c102^2 c021^2+36 c210^2 c111 c102 c021 c012+24 c210^2 c102^2 c030 c012-12 c210 c201^2 c030 c021 c003+24 c210 c201^2 c030 c012^2-12 c210 c201^2 c021^2 c012+6 c210 c201 c120^2 c012 c003-60 c210 c201 c120 c111 c021 c003+36 c210 c201 c120 c111 c012^2+18 c210 c201 c120 c102 c030 c003-6 c210 c201 c120 c102 c021 c012+36 c210 c201 c111^2 c030 c003-12 c210 c201 c111^2 c021 c012-60 c210 c201 c111 c102 c030 c012+36 c210 c201 c111 c102 c021^2+6 c210 c201 c102^2 c030 c021+12 c210 c120^2 c111 c102 c003-12 c210 c120^2 c102^2 c012-12 c210 c120 c111^3 c003-12 c210 c120 c111^2 c102 c012+36 c210 c120 c111 c102^2 c021-12 c210 c120 c102^3 c030+24 c210 c111^4 c012-36 c210 c111^3 c102 c021+12 c210 c111^2 c102^2 c030+4 c201^3 c030^2 c003-12 c201^3 c030 c021 c012+8 c201^3 c021^3+24 c201^2 c120^2 c021 c003-27 c201^2 c120^2 c012^2-24 c201^2 c120 c111 c030 c003+36 c201^2 c120 c111 c021 c012+6 c201^2 c120 c102 c030 c012-12 c201^2 c120 c102 c021^2+12 c201^2 c111^2 c030 c012-24 c201^2 c111^2 c021^2+12 c201^2 c111 c102 c030 c021-3 c201^2 c102^2 c030^2-12 c201 c120^3 c102 c003+12 c201 c120^2 c111^2 c003+36 c201 c120^2 c111 c102 c012-12 c201 c120^2 c102^2 c021-36 c201 c120 c111^3 c012-12 c201 c120 c111^2 c102 c021+12 c201 c120 c111 c102^2 c030+24 c201 c111^4 c021-12 c201 c111^3 c102 c030+8 c120^3 c102^3-8 c111^6;*)

(* The discriminant of the Hesse pencil is zero if and only if the cubic a*f + 6*hessian(f) degenerates into three lines. It is the cube of the following polynomial.*)

d = a^4 + 24 S a^2 + 8T a - 48 S^2;
sol = Solve[d == 0, a]; 

(* Let g be one of these degenerate cubics.  We factor g into three lines. *)

g = a f+6hessian/.sol[[3]];
monomials={x^2y,x y^2, y^3, x^2 z, xz^2, z^3, x  y  z, y^2 z, y z^2};
{A1,B1,C1,A2,B2,C2,D1, E1, F1}=Table[Coefficient[g,monomials[[i]]]/Coefficient[g,x^3],{i,1,9}];
generalcubic= Solve[z^3 + q2 z^2 + q1 z + q0 == 0,z];
asolutions = generalcubic/.{q0->-C1, q1->B1, q2-> -A1};

{a1,a2,a3}=z/.asolutions;
{b1,b2,b3}=(Inverse[{{s1,s2,s3},{t1,t2,t3},{u1,u2,u3}}]/.{s1->1,s2->1,s3->1,t1->a1,t2->a2,t3->a3,u1->a2 a3,u2->a1 a3,u3->a1 a2}).{A2, A1 A2-D1,E1};

(* Let h be another of these degenerate cubics.  We factor h into three lines. *)

h = a f+6hessian/.sol[[4]];
{A1,B1,C1,A2,B2,C2,D1, E1, F1}=Table[Coefficient[h,monomials[[i]]]/Coefficient[h,x^3],{i,1,9}];
asolutions = generalcubic/.{q0->-C1, q1->B1, q2-> -A1};
{c1,c2,c3}=z/.asolutions;
{d1,d2,d3}=(Inverse[{{s1,s2,s3},{t1,t2,t3},{u1,u2,u3}}]/.{s1->1,s2->1,s3->1,t1->c1,t2->c2,t3->c3,u1->c2 c3,u2->c1 c3,u3->c1 c2}).{A2, A1 A2-D1,E1};

(* Intersect the two sets of three lines pairwise to produce the nine inflection points. *)

inf = {
   {a1 d1-b1 c1, b1-d1, c1-a1},
   {a2 d1-b2 c1, b2-d1, c1-a2},
   {a3 d1-b3 c1, b3-d1, c1-a3},
   {a1 d2-b1 c2, b1-d2, c2-a1},
   {a2 d2-b2 c2, b2-d2, c2-a2},
   {a3 d2-b3 c2, b3-d2, c2-a3},
   {a1 d3-b1 c3, b1-d3, c3-a1},
   {a2 d3-b2 c3, b2-d3, c3-a2},
   {a3 d3-b3 c3, b3-d3, c3-a3}
   };
Table[{f,hessian}/.{x->inf[[i]][[1]],y->inf[[i]][[2]],z->inf[[i]][[3]]},{i,1,9}]

(* Write the result to a file. *)

inf>>"inf.txt"