/* InvPlaneQuartic.m

   Invariant Theory of plane quartic curves. 

   Shioda (1967) expects that the ring of invariants is generated by elements of degree at most 27. 
   We compute a list of invariants with the same degrees as given by Shioda.
   The completeness of this system is a proven fact.
 
   The algorithm requires characteristic > 5.

   Intrinsics in this package: 
       ContravariantsOfPlaneQuartic,  DixmierOhnoInv, DualCurve

   DualCurve is a general intrinsic to compute the (formal) dual curve of any plane curve.
   It illustrates the use of ClebschTransfer in a simple example.
 
*/

IsPlaneQuartic := function(f)
 if not IsHomogeneous(f) then return false; end if;
 if not Degree(f) eq 4 then return false; end if;
 if not Rank(Parent(f)) eq 3 then return false; end if;
 
 if Characteristic(Parent(f)) in [2,3,5] then return false; end if;

 return true;
end function;

/* Basic invariants of binary quartic. One of degree 2 and one of degree 3 
   The symbolic representations are (12)^4  and (12)^2 (13)^2 (23)^2 */
InvariantsOfBinaryQuartic := function(f, deg)
 r0 := Parent(f);
 a := MonomialCoefficient(f,r0.1^4);
 b := MonomialCoefficient(f,r0.1^3*r0.2);
 c := MonomialCoefficient(f,r0.1^2*r0.2^2); 
 d := MonomialCoefficient(f,r0.1^1*r0.2^3);
 e := MonomialCoefficient(f,       r0.2^4);
 if deg eq 2 then
  return 12*a*e - 3*b*d + c^2;
 end if;
 
 if deg eq 3 then
  return 72*a*c*e - 27*a*d^2 - 27*b^2*e + 9*b*c*d - 2*c^3;
 end if; 

 assert false;
end function;
 
intrinsic ContravariantsOfPlaneQuartic(f :: RngMPolElt) -> RngMPolElt,RngMPolElt
{The 2 contravariants of degree 2 and 3 (class 4 and 6) of a plane quartic curve that are 
 related to the basic invariants of a binary quartic. }

 require IsPlaneQuartic(f) : "Plane quartic in characteristic > 5 expected.";

 r0 := PolynomialRing(FieldOfFractions(BaseRing(Parent(f))),3);
 mon4 := Monomials((r0.1 + r0.2 + r0.3)^4);
 mon6 := Monomials((r0.1 + r0.2 + r0.3)^6);

 knots4 := [[1,i,j] : i,j in [0..4] | i+j le 4];

 knots6 := [[1,i,j] : i,j in [0..6] | i+j le 6];
 
 vcl2 := [NumericClebschTransfer(f, func<arg |  InvariantsOfBinaryQuartic(arg,2) >, knot) : knot in knots4];
 
 vcl3 := [NumericClebschTransfer(f, func<arg |  InvariantsOfBinaryQuartic(arg,3) >, knot) : knot in knots6];

 co2,k1 := Solution(Matrix([[Evaluate(mm,a) : a in knots4] : mm in mon4]), Vector(vcl2));
 co3,k2 := Solution(Matrix([[Evaluate(mm,a) : a in knots6] : mm in mon6]), Vector(vcl3));
 assert Dimension(k1) eq 0;
 assert Dimension(k2) eq 0;

 return &+[co2[i] * mon4[i] : i in [1..#mon4]], &+[co3[i] * mon6[i] : i in [1..#mon6]];
end intrinsic;

/* One complete system of invariants of plane quartic curves. */
intrinsic InvariantsOfPlaneQuartic(f :: RngMPolElt) -> SeqEnum
{A complete system of invariants of a plane quartic curve and the list of their degrees.}

 c_2_0_4,c_3_0_6 := ContravariantsOfPlaneQuartic(f);

 /* Hessian covaiant */
 c_3_6 := Determinant(Matrix([[Derivative(Derivative(f,i),j) : i in [1..3]] : j in [1..3]]));

 I_3 := ApplyContravariant(c_2_0_4,f);

 c_4_4 := ApplyContravariant(c_2_0_4,f^2); /* Qartic Covariant */
 I_6 := ApplyContravariant(c_2_0_4,c_4_4);

 c_4_0_2 := ApplyContravariant(f,c_3_0_6); /* Quadratic Contravariant */
 c_5_2 := ApplyContravariant(c_4_0_2,f);   /* Quadratic Covariant */
 I_9a := ApplyContravariant(c_4_0_2,c_5_2);

 c_7_4 := ApplyContravariant(c_2_0_4,c_4_4*f); /* Quartic Covariant */
 I_9b := ApplyContravariant(c_2_0_4,c_7_4);

 c_8_0_4 := ApplyContravariant(c_5_2,c_3_0_6); /* Quartic Contravariant */
 I_12a := ApplyContravariant(c_8_0_4, c_4_4);

 c_7_0_2 := ApplyContravariant(c_5_2,c_2_0_4);
 I_12b :=  ApplyContravariant(c_7_0_2, c_5_2);

 I_15a := ApplyContravariant(c_8_0_4, c_7_4);
 c_11_2 := ApplyContravariant(c_8_0_4, c_3_6);
 I_15b := ApplyContravariant(c_4_0_2, c_11_2);

 I_18a := ApplyContravariant(c_7_0_2,c_11_2); 
 c_13_0_2 := ApplyContravariant(c_5_2, c_8_0_4);
 I_18b := ApplyContravariant(c_13_0_2, c_5_2);
 
 c_14_2 := ApplyContravariant(c_7_0_2, c_7_4);
 I_21a := ApplyContravariant(c_7_0_2,c_14_2);
 c_17_2 := ApplyContravariant(c_13_0_2, c_4_4);
 I_21b := ApplyContravariant(c_4_0_2, c_17_2);

 c_20_2 := ApplyContravariant(c_13_0_2, c_7_4);
 I_27 := ApplyContravariant(c_7_0_2, c_20_2);

 return [MonomialCoefficient(a,1) : a in 
         [I_3,I_6,I_9a,I_9b, I_12a, I_12b, I_15a, I_15b, I_18a, I_18b, I_21a, I_21b, I_27]], 
        [3,6,9, 9, 12, 12, 15, 15, 18, 18, 21, 21, 27];
end intrinsic;

intrinsic DixmierOhnoInv(InvL :: SeqEnum) -> SeqEnum 
{Computes the Dixmier-Ohne invariants of a plane quartic.
 The 2nd D.-O. invariant is the catalectican. The last one is the discriminant.
 The input are the invariants as computed by InvariantsOfPlaneQuartic.
 [ 1, 2, 3, 5, 7, 9, 13] are the primary invariants of Dixmier.}

 I_3 := InvL[1];    I_6 := InvL[2];    I_9a := InvL[3];
 I_9b := InvL[4];   I_12a := InvL[5];  I_12b := InvL[6];
 I_15a := InvL[7];  I_15b := InvL[8];  I_18a := InvL[9];
 I_18b := InvL[10]; I_21a := InvL[11]; I_21b := InvL[12];
 I_27 := InvL[13];

return [1/6*I_3,-37/36*I_3^2 + 1/24*I_6, 1/3456*I_9a,
1295/1944*I_3^3 - 325/5184*I_3*I_6 + 1/31104*I_9a + 5/3456*I_9b,
1295/5184*I_3^4 - 325/13824*I_3^2*I_6 - 11/539136*I_3*I_9a + 5/9216*I_3*I_9b - 1/359424*I_12a + 1/22464
*I_12b,
-205/404352*I_3*I_9a + 5/269568*I_12a + 11/269568*I_12b,
2826985/46656*I_3^5 - 338365/41472*I_3^3*I_6 - 5605763/4852224*I_3^2*I_9a + 10915/82944*I_3^2*I_9b 
+ 19175/82944*I_3*I_6^2 + 33905/479232*I_3*I_12a + 23/89856*I_3*I_12b + 28787/995328*I_6*I_9a 
- 295/55296*I_6*I_9b - 79/36864 *I_15a + 1/55296*I_15b,
1747991/139968*I_3^5 - 210115/124416*I_3^3*I_6 - 795983/4852224*I_3^2*I_9a + 6781/248832*I_3^2*I_9b 
+ 11993/248832*I_3*I_6^2 + 383237/38817792*I_3*I_12a + 55/2426112*I_3*I_12b + 12349/2985984*I_6*I_9a 
- 185/165888*I_6*I_9b - 11/36864*I_15a + 1/497664*I_15b,
3951270589/929037600*I_3^6 + 6580481/20645280*I_3^4*I_6 - 29380644767/96619910400*I_3^3*I_9a 
- 2789051/103226400*I_3^3*I_9b - 449709877/6606489600*I_3^2*I_6^2 + 29011049/1590451200*I_3^2*I_12a 
+ 170771/1789257600*I_3^2*I_12b + 792900131/103061237760*I_3*I_6*I_9a + 2204177/440432640*I_3*I_6*I_9b 
- 8069/14681088*I_3*I_15a + 743/220216320*I_3*I_15b + 1039/146810880*I_6^3 - 331/159045120*I_6*I_12a 
- 2599/1272360960*I_6*I_12b + 335/264259584*I_9a^2 - 13/117448704*I_9a*I_9b - 1039/13049856*I_9b^2 
+ 1/7340544*I_18a - 67/264259584*I_18b,
673485492533/150504091200*I_3^6 - 15109137319/26756282880*I_3^4*I_6 
- 1170516456349/15652425484800*I_3^3*I_9a +  2181703373/267562828800*I_3^3*I_9b 
+ 7174470803/535125657600*I_3^2*I_6^2 + 20880416879/4637755699200*I_3^2*I_12a 
- 593863/289859731200*I_3^2*I_12b + 3921859229/2086990064640*I_3*I_6*I_9a 
- 2793689/17837521920*I_3*I_6*I_9b - 645371/4756672512*I_3*I_15a + 7849/8918760960*I_3*I_15b 
+ 21529/11891681280*I_6^3 - 39671/77295928320*I_6*I_12a + 158033/309183713280*I_6*I_12b 
+ 485/21405026304*I_9a^2 + 3257/9513345024*I_9a*I_9b - 3829/1057038336*I_9b^2 + 
    19/2378336256*I_18a - 97/21405026304*I_18b,
-38625671306177921/132443600256000*I_3^7 + 726062421254443/23545528934400*I_3^5*I_6 + 
    32677304297835107/5008776155136000*I_3^4*I_9a - 110610751507151/235455289344000*I_3^4*I_9b - 
    65180659422611/470910578688000*I_3^3*I_6^2 - 270649307084383/680204169216000*I_3^3*I_12a + 
    2221069283773/2040612507648000*I_3^3*I_12b - 35665690799237/459137814220800*I_3^2*I_6*I_9a - 
    514058092063/62788077158400*I_3^2*I_6*I_9b + 457742509/38053380096*I_3^2*I_15a 
- 2944371937/31394038579200*I_3^2*I_15b - 
    352578101171/20929359052800*I_3*I_6^3 - 5655504235193/1088326670745600*I_3*I_6*I_12a - 
    29464008091/272081667686400*I_3*I_6*I_12b - 44448449/7420409118720*I_3*I_9a^2 
- 88192739/6802041692160*I_3*I_9a*I_9b + 167337199/465096867840*I_3*I_9b^2 + 21911/38053380096*I_3*I_18a 
+ 7424489/3767284629504*I_3*I_18b - 609782189/283787919360*I_6^2*I_9a + 1222799/3153199104*I_6^2*I_9b 
+ 3365/21233664*I_6*I_15a - 60131/47297986560*I_6*I_15b + 893/10808328960*I_9a*I_12a 
+ 123241/691733053440*I_9a*I_12b - 8201/13974405120*I_9b*I_12a - 390833/153718456320*I_9b*I_12b +
    61/58392576*I_21a - 5/98537472*I_21b,
-247155462438880667/297998100576000*I_3^7 + 4417017123593911/52977440102400*I_3^5*I_6 + 
    53704049577589241/2817436587264000*I_3^4*I_9a - 693650140524377/529774401024000*I_3^4*I_9b + 
    10726661898583/66221800128000*I_3^3*I_6^2 - 10682973192629321/9182756284416000*I_3^3*I_12a + 
    527173712653/143480566944000*I_3^3*I_12b - 2971831691985509/16528961311948800*I_3^2*I_6*I_9a - 
    155267756693/4414786675200*I_3^2*I_6*I_9b + 6023560355/171240210432*I_3^2*I_15a 
- 19501945699/70636586803200*I_3^2*I_15b - 
    1458395910271/23545528934400*I_3*I_6^3 - 22051787886593/1224367504588800*I_3*I_6*I_12a - 
    3411484013/9565371129600*I_3*I_6*I_12b - 237988129/12521940387840*I_3*I_9a^2 
- 1591905703/30609187614720*I_3*I_9a*I_9b + 282680987/261616988160*I_3*I_9b^2 
+ 37577/21405026304*I_3*I_18a + 13130405/2119097604096*I_3*I_18b - 
    793041101/106420469760*I_6^2*I_9a + 1687493/1182449664*I_6^2*I_9b + 117853/214990848*I_6*I_15a 
- 76357/17736744960*I_6*I_15b + 48179/172933263360*I_9a*I_12a + 101917/172933263360*I_9a*I_12b 
- 2449/1746800640*I_9b*I_12a - 139897/19214807040*I_9b*I_12b + 157/49268736*I_21a - 1/6158592*I_21b,
-2349803391701124532313/3502855808075096417894400*I_3^9 
+ 357748057892475176519/15568248035889317412864000*I_3^7*I_6 + 
    676298725460507871857/44157212610886063934668800*I_3^6*I_9a 
- 1388242763499517637/2075766404785242321715200*I_3^6*I_9b + 
    20892893779275565961/4151532809570484643430400*I_3^5*I_6^2 
- 27997102082516662921/29983292513564611313664000*I_3^5*I_12a + 
    72686591753193881/29983292513564611313664000*I_3^5*I_12b + 
1332571116227484690043/1619097795732489010937856000*I_3^4*I_6*I_9a - 2169655194291454543/13838442698568282144768000*I_3^4*I_6*I_9b 
+ 78979914554477/2795644989609753968640*I_3^4*I_15a - 
    20251573641257/92256284657121880965120*I_3^4*I_15b 
- 950279304944345393/10378832023926211608576000* I_3^3*I_6^3 - 
    535368805171687003/7268676972985360318464000*I_3^3*I_6*I_12a 
- 16624676110057/169397132844997804032000*I_3^3*I_6*I_12b + 
    350574109536083/11480875278830376623013888*I_3^3*I_9a^2 
- 434843398325001823/53969926524416300364595200*I_3^3*I_9a*I_9b + 
    1489852697255671/922562846571218809651200*I_3^3*I_9b^2 
+ 1090660547/559128997921950793728*I_3^3*I_18a + 
    251710978865/55353770794273128579072*I_3^3*I_18b 
- 7895831234302641781/431759412195330402916761600*I_3^2*I_6^2*I_9a + 
    2059928851410157/922562846571218809651200*I_3^2*I_6^2*I_9b 
+ 1107969866929/497003553708400705536*I_3^2*I_6*I_15a - 
    3644780639399/205013965904715291033600*I_3^2*I_6*I_15b 
- 4986950714574397/519710403568453262770176000*I_3^2*I_9a*I_12a + 
    649714873372039/64963800446056657846272000*I_3^2*I_9a*I_12b +
 197946483914657/403815387388075573248000*I_3^2*I_9b*I_12a - 
    1396478329/135558143959620096000*I_3^2*I_9b*I_12b + 29406433/7721806625413005312*I_3^2*I_21a - 
    33013/285992837978259456*I_3^2*I_21b - 55580234364107/22779329544968365670400*I_3*I_6^4 - 
    12811055552742631/17767877045075325222912000*I_3*I_6^2*I_12a 
- 70821762769321/4441969261268831305728000*I_3*I_6^2*I_12b - 
    3340674864079/2355051339247256743182336*I_3*I_6*I_9a^2 
+ 2833678683412529/14391980406511013430558720*I_3*I_6*I_9a*I_9b + 
    2611260647861/123008379542829174620160*I_3*I_6*I_9b^2 + 2343317/41416962809033392128*I_3*I_6*I_18a + 
    1768954481/12300837954282917462016*I_3*I_6*I_18b + 4333377359/48457846486569068789760*I_3*I_9a*I_15a - 
    70824449/36343384864926801592320*I_3*I_9a*I_15b - 6144321929/414169628090333921280*I_3*I_9b*I_15a +
    18920833/155313610533875220480*I_3*I_9b*I_15b + 3963276137/9886252421930287104000*I_3*I_12a^2 - 
    8697667139/14829378632895430656000*I_3*I_12a*I_12b - 35974997/7414689316447715328000*I_3*I_12b^2 - 
    48805034529727/164011172723772232826880*I_6^3*I_9a + 512489707913/9111731817987346268160*I_6^3*I_9b
 +  615415243/28079296819683655680*I_6^2*I_15a - 8966651/51478710836086702080*I_6^2*I_15b + 
    147724642771/888393852253766261145600*I_6*I_9a*I_12a 
- 109101767183/444196926126883130572800*I_6*I_9a*I_12b - 
    188206219/4486837637645284147200*I_6*I_9b*I_12a - 4774514911/49355214014098125619200*I_6*I_9b*I_12b
 +  428633/5147871083608670208*I_6*I_21a - 377/190661891985506304*I_6*I_21b 
- 979/6368828741954720759808*I_9a^3 + 
    14817821/30192965887785342861312*I_9a^2*I_9b + 133705799/6709547975063409524736*I_9a*I_9b^2 - 
    1/5824260395020320768*I_9a*I_18a + 1/26209171777591443456*I_9a*I_18b 
- 521660159/2236515991687803174912*I_9b^3 - 
    3581/41416962809033392128*I_9b*I_18a - 31363/46594083160162566144*I_9b*I_18b 
- 362887/30419238221323960320*I_12a*I_15a + 
    18929/228144286659929702400*I_12a*I_15b + 264707/15209619110661980160*I_12b*I_15a 
- 6751/57036071664982425600*I_12b*I_15b +  1/2888816545234944*I_27];
end intrinsic;


DiscriminantOfBinaryForm := function(f)
 deg := Degree(f); 
 r := Parent(f);
 coeff := [ MonomialCoefficient(f,r.2^i * r.1^(deg - i)) : i in [0..deg]];
 
 di := Derivative(f,1);
 coeff_di :=  [ MonomialCoefficient(di,r.2^i * r.1^(deg - i - 1)) : i in [0..deg-1]];
 
 mat := Matrix([ [0 : tt in [2..i]] cat coeff cat [0 : tt in [i..deg-2]] : i in [1..deg-1]] cat 
               [ [0 : tt in [2..i]] cat coeff_di cat [0 : tt in [i..deg-1]] : i in [1..deg]]);
 mat[1,1] := 1; mat[deg,1] := deg;
 return Determinant(mat);
end function;

intrinsic DualCurve(f :: RngMPolElt) -> RngMPolElt 
{Computes the dual curve of the curve f as a contravariant of degree deg(f)(deg(f) - 1).}
 
 require (IsHomogeneous(f) and (Rank(Parent(f)) eq 3) and (Degree(f) gt 1)) : 
           "Curve must be given as a homogeneous polynomial in 3 variables of degree at least 2.";

 deg := Degree(f);
 deg_dual := deg*(deg - 1); 

 require ((deg_dual lt Characteristic(Parent(f))) or (Characteristic(Parent(f)) eq 0)) 
         : "Characteristic must be bigger than the degree of the result.";

 r := PolynomialRing(FieldOfFractions(BaseRing(Parent(f))),3);
 mon := [r.1^e1 * r.2^e2 * r.3^(deg_dual - e1 - e2) : e1,e2 in [0..deg_dual] | e1 +e2 le deg_dual];
 
 knots := [[1,j,k] : j,k in [0..deg_dual] | j+k le deg_dual+1];

 mat :=  Matrix([[Evaluate(m,pt) : pt in knots] : m in mon]);
 vec := Vector([ NumericClebschTransfer(f, DiscriminantOfBinaryForm, pt) : pt in knots]);
// printf"Solving linear system\n";
 sol,ker := Solution(mat,vec);
 assert Dimension(ker) eq 0;

 return &+[sol[i] * mon[i] : i in [1..#mon]]; 
end intrinsic;

