// ***************************************************************
// ****** This is the file precomputation_5.m accompanying *******
// *** LIFTING LOW-GONAL CURVES FOR USE IN TUITMAN'S ALGORITHM ***
// *********** by Wouter Castryck and Floris Vermeulen ***********
// ***************************************************************

RR<a,b,c,d,e> := PolynomialRing(Rationals(),5);
R<y1,y2,y3,y4,y5> := PolynomialRing(FieldOfFractions(RR),5, "grevlexw", [1,50,50,50,50]);
Q, pi := quo< R | y1^5 + a*y1^4 + b*y1^3 + c*y1^2 + d*y1 + e,
                  y2^5 + a*y2^4 + b*y2^3 + c*y2^2 + d*y2 + e,
                  y3^5 + a*y3^4 + b*y3^3 + c*y3^2 + d*y3 + e,
                  y4^5 + a*y4^4 + b*y4^3 + c*y4^2 + d*y4 + e,
                  y5^5 + a*y5^4 + b*y5^3 + c*y5^2 + d*y5 + e,
                  ElementarySymmetricPolynomial(R, 1) + a,
                  ElementarySymmetricPolynomial(R, 2) - b,
                  ElementarySymmetricPolynomial(R, 3) + c,
                  ElementarySymmetricPolynomial(R, 4) - d,
                  ElementarySymmetricPolynomial(R, 5) + e>;

sqrtdiscf := Determinant(Matrix(R, 5, 5,
             [    1,    1,    1,    1,    1,
                 y1,   y2,   y3,   y4,   y5,
               y1^2, y2^2, y3^2, y4^2, y5^2,
               y1^3, y2^3, y3^3, y4^3, y5^3,
               y1^4, y2^4, y3^4, y4^4, y5^4 ]));

rhos := [(y1*y2 + y2*y3 + y3*y4 + y4*y5 + y5*y1 - y1*y3 - y3*y5 - y5*y2 - y2*y4 - y4*y1)^2,
         (y1*y2 + y2*y4 + y4*y3 + y3*y5 + y5*y1 - y1*y3 - y3*y2 - y2*y5 - y5*y4 - y4*y1)^2,
         (y1*y2 + y2*y3 + y3*y5 + y5*y4 + y4*y1 - y1*y3 - y3*y4 - y4*y2 - y2*y5 - y5*y1)^2,
         (y1*y2 + y2*y4 + y4*y5 + y5*y3 + y3*y1 - y1*y4 - y4*y3 - y3*y2 - y2*y5 - y5*y1)^2,
         (y1*y2 + y2*y5 + y5*y4 + y4*y3 + y3*y1 - y1*y4 - y4*y2 - y2*y3 - y3*y5 - y5*y1)^2,
         (y1*y2 + y2*y5 + y5*y3 + y3*y4 + y4*y1 - y1*y3 - y3*y2 - y2*y4 - y4*y5 - y5*y1)^2];

rhos := [pi(rho) : rho in rhos];

// write hard-coded function for computing Cayley's sextic resolvent

print "Computing Cayley's sextic resolvent...";
tim := Cputime();

Rres<rho1,rho2,rho3,rho4,rho5,rho6> := PolynomialRing(Q,6);
S<x> := PolynomialRing(RR);

resolvent := x^6
             - (RR ! Evaluate(ElementarySymmetricPolynomial(Rres,1),rhos))*x^5
             + (RR ! Evaluate(ElementarySymmetricPolynomial(Rres,2),rhos))*x^4
             - (RR ! Evaluate(ElementarySymmetricPolynomial(Rres,3),rhos))*x^3
             + (RR ! Evaluate(ElementarySymmetricPolynomial(Rres,4),rhos))*x^2
             - (RR ! Evaluate(ElementarySymmetricPolynomial(Rres,5),rhos))*x
             + (RR ! Evaluate(ElementarySymmetricPolynomial(Rres,6),rhos));

Write("precomputed_5.m", "CayleyResolvent := function(monic_pol);" : Overwrite := true);
Write("precomputed_5.m", "  R<x> := Parent(monic_pol);");
Write("precomputed_5.m", "  a := Coefficient(monic_pol,4);");
Write("precomputed_5.m", "  b := Coefficient(monic_pol,3);");
Write("precomputed_5.m", "  c := Coefficient(monic_pol,2);");
Write("precomputed_5.m", "  d := Coefficient(monic_pol,1);");
Write("precomputed_5.m", "  e := Coefficient(monic_pol,0);");
Write("precomputed_5.m", "  return ");
Write("precomputed_5.m", resolvent);
Write("precomputed_5.m", ";");
Write("precomputed_5.m", "end function;\n");

print "This took", Cputime(tim);

// write hard-coded function for computing altered fundamental resolvent map

V, vecpi := VectorSpace(Q);
M := [ZeroMatrix(RR, 5, 5) : i in [0..4]];

print "Computing the fundamental resolvent map...";
tim := Cputime();

for i in [1..5] do
//  for r in [1..5] do
//    M[r][i,i] := 0;
//  end for;
  for j in [i+1..5] do
    print "Doing entry", i, j, "...";
    entry := pi(sqrtdiscf) * Determinant(Matrix(Q, 3, 3,
             [                     1,                     1,                     1,
               rhos[1]^i + rhos[2]^i, rhos[3]^i + rhos[6]^i, rhos[4]^i + rhos[5]^i,
               rhos[1]^j + rhos[2]^j, rhos[3]^j + rhos[6]^j, rhos[4]^j + rhos[5]^j ]));
    entry := vecpi(entry);
    for r in [1..5] do
      M[r][i,j] := RR ! entry[r];
//      M[r][j,i] := -M[r][i,j];
    end for;
  end for;
end for;

Write("precomputed_5.m", "FundamentalResolventMatrices := function(monic_pol);");
Write("precomputed_5.m", "  R<x> := Parent(monic_pol);");
Write("precomputed_5.m", "  a := Coefficient(monic_pol,4);");
Write("precomputed_5.m", "  b := Coefficient(monic_pol,3);");
Write("precomputed_5.m", "  c := Coefficient(monic_pol,2);");
Write("precomputed_5.m", "  d := Coefficient(monic_pol,1);");
Write("precomputed_5.m", "  e := Coefficient(monic_pol,0);");
Write("precomputed_5.m", "  M := [ZeroMatrix(BaseRing(R),5,5) : i in [0..4]];");
for i in [1..5] do
  for j in [i+1..5] do
    for r in [1..5] do
      Write("precomputed_5.m", "M["*IntegerToString(r)*"]["*IntegerToString(i)*","*IntegerToString(j)*"] :=");
      Write("precomputed_5.m", M[r][i][j]);
      Write("precomputed_5.m", ";");
      Write("precomputed_5.m", "M["*IntegerToString(r)*"]["*IntegerToString(j)*","*IntegerToString(i)*"] :=
                                -M["*IntegerToString(r)*"]["*IntegerToString(i)*","*IntegerToString(j)*"];");
    end for;
  end for;
end for;
Write("precomputed_5.m", "  return M;");
Write("precomputed_5.m", "end function;");
print "This took", Cputime(tim);
