// THIS FILE ACCOMPANIES THE PAPER "Scrollar invariants, syzygies and representations of
// the symmetric group" BY W. CASTRYCK, F. VERMEULEN AND Y. ZHAO. IT COMPUTES VALUES OF
// THE SCROLLAR INVARIANTS OF RESOLVENT CURVES OF DEGREE 6 COVERS, USING THE INTRINSIC
// ShortBasis(). THE OUTPUT LENDS SUPPORT FOR CONJECTURE 44.

// WARNING: THIS CODE NEEDS FURTHER DOCUMENTATION AND CLEANING UP (IN PROGRESS)!

// INPUT. DEGS DESCRIBES NEWTON POLYGON, SHOULD BE CONVEX OF LENGTH 6.
// THE GROUP INDEX CORRESPONDS TO THE SAGEMATH INDEX FROM OUR PAPER.

aboutprimesize := 2^16;
degs := [13,11,9,7,5,3];
t_adic_precision := 2000;
group_index := 33;

// ********************************************************
// GENERATE f, COMPUTE ITS MARONI AND SCHREYER INVARIANTS *
// ********************************************************

load "maroni_inv.m";

p := NextPrime(aboutprimesize);
F := GF(p);

R<t> := PolynomialRing(F);
RR<T> := PowerSeriesRing(F,t_adic_precision);
complete := hom<R -> RR | T>;
chop := hom<RR -> R | t>;

S<x> := PolynomialRing(R);
SS<X> := PolynomialRing(RR);

f := &*[x - a : a in [1..6]]; // DISTINCT ROOTS MOD t
for i in [0..5] do
  if degs[i+1] ne 0 then
    f +:= &+[Random(F)*t^j : j in [1..degs[i+1]]]*x^i;
  end if;
end for;
coeffs_f := Coefficients(f);

// COMPUTE ITS MARONI AND SCHREYER INVARIANTS

bivarR := PolynomialRing(F,2);
emb := hom<R -> bivarR | bivarR.1>;
f_bivar := &+[emb(coeffs_f[i])*bivarR.2^(i-1) : i in [1..#coeffs_f]];

Delta := NewtonPolytope(f_bivar);
Delta1 := InnerPoints(Delta);

g := #Delta1;

print "Computed a randomish polynomial f with Newton polygon";
print Delta; 
print "  (it should define a curve of genus", g, ", which is",
      g eq Genus(Curve(AffinePlane(F),f_bivar)), "-- if false then everything is meaningless).";

e_i := MaroniInvariants(f_bivar,1);
// print "Its Maroni invariants are", e_i;
b_i := [b + 2 : b in SchreyerInvariants(Delta)];
// print "while its augmented Schreyer invariants are", b_i;

// COMPUTE ITS ROOTS AND AN EVALUATION MAP FROM F[alp1, ..., alp6]

fcomplete := &+[complete(coeffs_f[i+1])*X^i : i in [0..6]];
rts := Roots(fcomplete); // SHOULD HAVE 6 DISTINCT ROOTS BY HENSEL'S LEMMA
Rund_rts<alp1, alp2, alp3, alp4, alp5, alp6> := PolynomialRing(RR,6);
evmap := hom<Rund_rts -> RR | [rts[i][1] : i in [1..6]] >;

print "Computed its approximated roots, now ready to compute resolvents...";

// ****************************
// CONSTRUCT EXOTIC RESOLVENT *
// ****************************

S6 := SymmetricGroup(6);
exoticgp := sub<S6 | S6 ! (1,2,3,4), S6 ! (1,5,6,2) >;
to_be_permuted := &*[Rund_rts.i^(i-1) : i in [1..6]];
inv_expression := Rund_rts ! 0;
for sigma in exoticgp do
  indices := Eltseq(sigma);
  inv_expression +:= Evaluate(to_be_permuted, [Rund_rts.indices[i] : i in [1..6]]);
end for;

conjugates := {inv_expression};
for pi in Permutations({1..6}) do
  conjugates join:= {Evaluate(inv_expression,[Rund_rts.i : i in pi])};
end for;
exotic_resolvent := &*[X - evmap(c) : c in conjugates];

// NOW CHOP COEFFICIENTS;
coeffs_ex_res := [chop(c) : c in Coefficients(exotic_resolvent)];
exotic_resolvent := &+[coeffs_ex_res[i+1]*x^i : i in [0..Degree(exotic_resolvent)]];

print "   ";
print "Computed its exotic resolvent, which should be of degree 6, which is", #conjugates eq 6;
exotic_resolvent_bivar := &+[emb(coeffs_ex_res[i])*bivarR.2^(i-1) : i in [1..#coeffs_ex_res]];
a_i := MaroniInvariants(exotic_resolvent_bivar,1);
// print "It has Maroni invariants (i.e., the invariants a_1, ..., a_5) are", a_i;

// ******************************
// COMPUTE RECURRING INVARIANTS *
// ******************************

e_i2 := [];
for i in [1..4] do for j in [i+1..5] do
  e_i2 cat:= [e_i[i] + e_i[j] + 2];
end for; end for;
e_i2 := Sort(e_i2);

e_i3 := [];
for i in [1..3] do for j in [i+1..4] do for k in [j+1..5] do
  e_i3 cat:= [e_i[i] + e_i[j] + e_i[k] + 4];
end for; end for; end for;
e_i3 := Sort(e_i3);

e_i4 := Sort([g+1-e : e in e_i]);

b_i3 := Sort([g+1-b : b in b_i]);

ga_i := Sort([g+1-a : a in a_i]);

print "   ";
print "The pool of invariants consists of:";
print " --> the e_i", e_i; 
print "     they should sum to", g - 5, "which is", &+e_i eq g - 5;
print " --> the e_i^(2) + 2", e_i2;
print "     they should sum to", 4*g, "which is", &+e_i2 eq 4*g;
print " --> the e_i^(3) + 4", e_i3;
print "     they should sum to", 6*g + 10, "which is", &+e_i3 eq 6*g + 10;
print " --> the e_i^(4) + 6", e_i4;
print "     they should sum to", 4*g + 10, "which is", &+e_i4 eq 4*g + 10;
print " --> g + 3 which is", g+3;
print " --> the b_i + 2", b_i;
print "     they should sum to", 3*g - 3, "which is", &+b_i eq 3*g - 3;
print " --> the b_i^(2) + 4, which we don't know explicitly";
print "     these are 16 numbers which pairwise sum up to", g + 1, "for a total of", 8*g+8;
print " --> the b_i^(3) + 6", b_i3;
print "     they should sum to", 6*g + 12, "which is", &+b_i3 eq 6*g + 12;
print " --> the exotic invariants a_i", a_i;
print "     they should sum to", 3*g + 5, "which is", &+a_i eq 3*g + 5;
print " --> the g + 1 - a_i", ga_i;
print "     they should sum to", 2*g, "which is", &+ga_i eq 2*g;

// TODO, MAYBE: INCORPORATE GUESSMULTS IN HERE

// ******************************
// CONSTRUCT SPECIFIC RESOLVENT *
// ******************************

print "   ";

case group_index:
  when 32:
    gp := sub<S6 | S6 ! (1,2), S6! (1,2,3), S6 ! (4,5)>;
    index := 60;
    genus_res := 22*g + 51;
    inv_expression := alp1 + alp2 + alp3 - alp4 - alp5;
    conjecturetext := "no conjecture yet";
    conjecture := "to do";
    completeconjecture := false;
  when 33:
    gp := sub<S6 | S6 ! (1,2), S6 ! (1,3,2,4), S6 ! (5,6) >;
    index := 45;
    genus_res := 18*g + 46;
    inv_expression := alp1*alp2 + alp3*alp4;
    conjecturetext := "e_i, b_i + 2 (2x), b_i^(2) + 4, a_i";
    conjecture := Sort(e_i cat b_i cat b_i cat a_i);
    completeconjecture := false;
  when 35:
    gp := sub<S6 | S6 ! (1,2,3,4,5,6), (S6 ! (1,5,3))*(S6 ! (4,6,2))>;
    index := 40;
    genus_res := 20*g + 61;
    inv_expression := Rund_rts ! 0;
    conjecturetext := "no conjecture yet, according to orbits a_i and b_i + 2, no e_i or g + 3";
    conjecture := "to do";
    completeconjecture := false;
  when 36:
    gp := sub<S6 | S6 ! (4,5), (S6 ! (1,2,3))*(S6 ! (4,5)), S6 ! (4,5,6), (S6 ! (1,2,3))*(S6 ! (4,5,6)) >;
    index := 40;
    genus_res := 16*g + 41;
    inv_expression := Rund_rts ! 0;
    conjecturetext := "no conjecture yet, based on orbits expect e_i, b_i+2, no a_i's and no g+3";
    conjecture := "to do";
    completeconjecture := false; 
  when 38:
    gp := sub<S6 | S6 ! (1,2), S6 ! (3,4,5), S6 ! (3,4,6), (S6 ! (3,4))*(S6 ! (5,6))>;
    index := 30;
    genus_res := 14*g + 41;
    inv_expression := Rund_rts ! 0;  
    conjecturetext := "e_i, e_i^(3) + 4, e_i^(4) + 6, b_i + 2";
    conjecture := Sort(e_i cat e_i3 cat e_i4 cat b_i);
    completeconjecture := true;
  when 40:
    gp := sub<S6 | (S6 ! (1,4,5))*(S6 ! (2,3,6)), S6 ! (1,6,2,5) >;
    index := 30;
    genus_res := 15*g + 46;
    inv_expression := Rund_rts ! 0;
    conjecturetext := "no conjecture yet, based on orbits expect a_i (2x), b_i + 2, no e_i and no g + 3";
    conjecture := "to do";
    completeconjecture := false;
  when 41:
    gp := sub<S6 | S6 ! (1,2), S6 ! (2,3), S6 ! (3,4), S6 ! (4,1) >;
    index := 30;
    genus_res := 9*g+16;
    inv_expression := alp5 - alp6;
    conjecturetext := "e_i (2x), e_i^(2) + 2, b_i + 2";
    conjecture := Sort(e_i cat e_i cat e_i2 cat b_i);
    completeconjecture := true;
  when 45:
    gp := sub<S6 | S6 ! (1,2), S6 ! (1,2,3), S6 ! (4,5), S6 ! (4,5,6) >;
    index := 20;
    genus_res := 6*g+11;
    inv_expression := alp1 + alp2 + alp3;
    conjecturetext := "e_i, b_i + 2, g + 1 - a_i";
    conjecture := Sort(e_i cat b_i cat ga_i);
    completeconjecture := true;
  when 49:
    gp := sub<S6 | S6 ! (1,2,3), S6 ! (3,4,5) >;
    index := 12;
    genus_res := 6*g + 19;
    inv_expression := (alp1 - alp2)*(alp1 - alp3)*(alp1 - alp4)*(alp1 - alp5)*
                      (alp2 - alp3)*(alp2 - alp4)*(alp2 - alp5)*
                      (alp3 - alp4)*(alp3 - alp5)*(alp4 - alp5);
    conjecturetext := "e_i, e_i^(4) + 6, g + 3";
    conjecture := Sort(e_i cat e_i4 cat [g+3]);
    completeconjecture := true;
  when 54:
    gp := AlternatingGroup(6);
    index := 2;
    genus_res := g + 4;
    inv_expression := (alp1 - alp2)*(alp1 - alp3)*(alp1 - alp4)*(alp1 - alp5)*(alp1 - alp6)*
                      (alp2 - alp3)*(alp2 - alp4)*(alp2 - alp5)*(alp2 - alp6)*
                      (alp3 - alp4)*(alp3 - alp5)*(alp3 - alp6)*
                      (alp4 - alp5)*(alp4 - alp6)*(alp5 - alp6);
    conjecturetext := "g+3";
    conjecture := [g+3];
    completeconjecture := true;
end case;

print "Now turning specifically to the resolvent for the group";
print gp;
print "We conjecture that the corresponding scrollar invariants are: "*conjecturetext;

if inv_expression eq 0 then
  // AUTOMATED SEARCH FOR AN INVARIANT EXPRESSION, TRY TO KEEP DEGREE LOW
  for numvars in [2..5] do
    for i in [1..10] do // TEN ATTEMPTS PER numvars
      pi := Random(Permutations({1..6}));
      to_be_permuted := &*[Rund_rts.pi[i]^i : i in [1..numvars]];
      inv_expression := Rund_rts ! 0;
      for sigma in gp do
        inv_expression +:= Evaluate(to_be_permuted, [Rund_rts.Eltseq(sigma)[i] : i in [1..6]]);
      end for;
      conjugates := {inv_expression};
      for pi in Permutations({1..6}) do
        conjugates join:= {Evaluate(inv_expression,[Rund_rts.i : i in pi])};
      end for;
      if #conjugates eq index then break numvars; end if;
    end for;
  end for;
  print "Found invariant expression of degree", Degree(to_be_permuted);
else
  conjugates := {inv_expression};
  for pi in Permutations({1..6}) do
    conjugates join:= {Evaluate(inv_expression,[Rund_rts.i : i in pi])};
  end for;
end if;

resolvent := &*[X - evmap(c) : c in conjugates];

// NOW CHOP COEFFICIENTS;
coeffs_res := [chop(c) : c in Coefficients(resolvent)];
resolvent := &+[coeffs_res[i+1]*x^i : i in [0..Degree(resolvent)]];
print "Computed the resolvent for this group, which should be of degree", index, "which is", #conjugates eq index;
print "Its expected genus is", genus_res;

// ******************************************************************
// FINALLY, COMPUTE THE MARONI INVARIANTS AND COMPARE TO CONJECTURE *
// ******************************************************************

resolvent_bivar := &+[emb(coeffs_res[i])*bivarR.2^(i-1) : i in [1..#coeffs_res]];
final_e_i := MaroniInvariants(resolvent_bivar,1);
print "Its Maroni invariants are", final_e_i;
if completeconjecture 
  then 
    print "Are these as conjectured?", final_e_i eq conjecture;
  else
    print "to do", final_e_i;
end if;
