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

load "precomputed_5.m"; // load precomputed Cayley resolvent and fundamental resolvent map

// PARAMETERS

F<a> := FiniteField(11,2);
R<x> := PolynomialRing(F);
g := 7;
e := [0,1,1,1]; // MARONI INVS, SUM SHOULD BE g - 4;
b := [0,0,0,1,1]; // SCHREYER INVS, SUM SHOULD BE 2*g - 12;

print "----";
print "Testing our lifting procedure for a generic five-gonal curve of genus", g;
print "with Maroni invariants", e[1], e[2], e[3], e[4], "and Schreyer invariants", b[1], b[2], b[3], b[4], b[5];
print "The Maroni invariants should sum up to g - 4, which is", &+e eq g-4;
print "The Schreyer invariants should sum up to 2g-12, which is", &+b eq 2*g-12;
print "----";

if Degree(F) eq 1 then
    K := Rationals();
else
    K<a> := ext<Rationals() | PolynomialRing(Integers()) ! MinimalPolynomial(a)>;
end if;

// AUXILIARY FUNCTIONS FOR LIFTING

NaiveLiftEltModp := function(alpha)
  p := #Parent(alpha);
  lift := Integers() ! alpha;
  if lift gt p/2 then return lift - p; else return lift; end if;
end function;

NaiveLiftCoeff := function(alpha, K)
  n := Degree(K);
  if n eq 1 then
    return K ! NaiveLiftEltModp(alpha);
  else
    elt := Eltseq(alpha);
    return &+[NaiveLiftEltModp(elt[i])*K.1^(i-1) : i in [1..n]];
  end if;
end function;

NaiveLiftPol := function(pol, K);
  d := Degree(pol);
  return PolynomialRing(K) ! [NaiveLiftCoeff(alpha, K) : alpha in Eltseq(pol)]; 
end function;

// GENERATE RANDOM DEGREE 5 COVER WITH ABOVE PARAMETERS
// AND SCRAMBLE OUTPUT SOMEWHAT

RandomSkewSymmetricMatrix := function(F, e, b)
  g := &+e + 4;
  R<x> := PolynomialRing(F);
  S<y1,y2,y3> := PolynomialRing(R,3);
  M := ZeroMatrix(S, 5, 5);
  for i in [1..5] do
    for j in [i+1..5] do
      if #[0..e[4] + b[i] + b[j] + 6 - g] ge 1 then M[i][j] +:= &+[Random(F)*x^ell : ell in [0..e[4] + b[i] + b[j] + 6 - g]]; end if;
      if #[0..e[1] + b[i] + b[j] + 6 - g] ge 1 then M[i][j] +:= &+[Random(F)*x^ell : ell in [0..e[1] + b[i] + b[j] + 6 - g]]*y1; end if;
      if #[0..e[2] + b[i] + b[j] + 6 - g] ge 1 then M[i][j] +:= &+[Random(F)*x^ell : ell in [0..e[2] + b[i] + b[j] + 6 - g]]*y2; end if;
      if #[0..e[3] + b[i] + b[j] + 6 - g] ge 1 then M[i][j] +:= &+[Random(F)*x^ell : ell in [0..e[3] + b[i] + b[j] + 6 - g]]*y3; end if;
      M[j][i] := -M[i][j];
    end for;
  end for;
  return M;
end function;

S<y> := PolynomialRing(R);
SS<y1,y2,y3> := PolynomialRing(R,3);
phi := hom<SS -> S | 0, 0, y>;

M := RandomSkewSymmetricMatrix(F, e, b);
f := phi(Basis(EliminationIdeal(Ideal(Pfaffians(M, 4)),2))[1]);
f := y^5 + Coefficient(f,4)*y^4
     + Coefficient(f,3)*Coefficient(f,5)*y^3
     + Coefficient(f,2)*Coefficient(f,5)^2*y^2
     + Coefficient(f,1)*Coefficient(f,5)^3*y
     + Coefficient(f,0)*Coefficient(f,5)^4; // MAKE MONIC
f := Evaluate(f, y + x^10 + x); // SOME RANDOM FURTHER SCRAMBLING

// CHECK IF RANDOM CURVE HAS THE CORRECT GENUS

FFf := FunctionField(f); gf := Genus(FFf);
print "Genus of input curve is", gf,"which equals prescribed genus g:", g eq gf;
print "----";

// COMPUTE AND LIFT THE SKEW-SYMMETRIC MATRIX

dualBasis := function(b)
  K := Parent(b[1]);
  frkx := BaseField(K);
  n := #b;
  C := ZeroMatrix(frkx, n^2, n^2);
  w := Vector(frkx, [0 : i in [1..n^2]]);
  for i in [1..n] do
    for j in [1..n] do
      for l in [1..n] do
        C[n*(i-1)+j, n*(j-1)+l] := Trace(b[i]*b[l]);
        if i eq j then w[n*(i-1)+j] := 1; end if; 
      end for;
    end for;
  end for;
  sol := Solution(C,w);
  db := [&+[sol[n*(j-1)+l]*b[l] : l in [1..n]] : j in [1..n]];
  return db;
end function;

BhargavaMatrices := function(f);
  kx := BaseRing(Parent(f)); frkx := FieldOfFractions(kx);
  fres := CayleyResolvent(f);
  FF<Y> := FunctionField(f); FFres<Yres> := FunctionField(fres);
  gfres := Genus(FFres);
  // if following test fails, cover is not simply branched
  print "Genus of resolvent curve is", gfres, "which equals 3g + 7:", gfres eq 3*gf+7;
  print "----";

  // computing altered fundamental resolvent map
  Mo := FundamentalResolventMatrices(f);
  print "Computed altered fundamental resolvent matrices";

  // computing reduced bases and their duals
  sb, hess := ShortBasis(DivisorGroup(FF) ! 0);
  sbres, hessres := ShortBasis(DivisorGroup(FFres) ! 0);
  dsb := dualBasis(sb);
  dsbres := dualBasis(sbres);

  b := [1,Y,Y^2,Y^3,Y^4]; //naive basis for FF
  N, h := Module(dsb, frkx : IsBasis := true);
  bres := [1,Yres,Yres^2,Yres^3,Yres^4,Yres^5]; //naive basis for FFres
  Nres, hres := Module(bres, frkx : IsBasis := true);

  //Computing the constant factor
  Ndisc, hdisc := Module(b, frkx : IsBasis := true);
  discmatr := Matrix(frkx, 5,5, [Inverse(hdisc)(sb[i])[j] : i in [1,2,3,4,5], j in [1,2,3,4,5]]); //change of basis from sb to b
  const := Determinant(discmatr); //The determinant of this matrix is (disc k[C]_0 / disc f)^(1/2)

  //Step 4: change of basis matrix
  A := Matrix(frkx, 5,5, [Inverse(h)(b[i])[j] : i in [1,2,3,4,5], j in [1,2,3,4,5]]); //Change of basis from b to dsb
  B := Matrix(frkx, 5,5, [Inverse(hres)(dsbres[i])[j] : i in [2,3,4,5,6], j in [2,3,4,5,6]]); //change of basis from dsbres to bres

  //Step 5: apply the action
  basechangedMo := [ZeroMatrix(frkx, 5,5) : i in [0..4]];
  for i in [1..5] do
    basechangedMo[i] := Transpose(B)*Mo[i]*B;
  end for;

  newMo := [ZeroMatrix(frkx, 5,5) : i in [0..4]];
  for i in [1..5] do
    for j in [1..5] do
      newMo[i] := newMo[i] + const*A[i,j]*basechangedMo[j];
    end for;
  end for;

  print "Applied base change to obtain genuine fundamental resolvent matrices";
  print "----";

  //The model is a Pfaffian model defined by newMo[2],newMo[3],newMo[4],newMo[5]
  kx5x5 := Parent(ZeroMatrix(kx,5,5));

  return kx5x5 ! newMo[2],
         kx5x5 ! newMo[3],
         kx5x5 ! newMo[4],
         kx5x5 ! newMo[5];
end function;

M1, M2, M3, M4 := BhargavaMatrices(f);

Rlift<X> := PolynomialRing(K);

M1lift := ZeroMatrix(Rlift, 5, 5);
for i in [1..5] do for j in [i+1..5] do
  M1lift[i][j] := NaiveLiftPol(M1[i][j], K);
  M1lift[j][i] := -M1lift[i][j];
end for; end for;
M2lift := ZeroMatrix(Rlift, 5, 5);
for i in [1..5] do for j in [i+1..5] do
  M2lift[i][j] := NaiveLiftPol(M2[i][j], K);
  M2lift[j][i] := -M2lift[i][j];
end for; end for;
M3lift := ZeroMatrix(Rlift, 5, 5);
for i in [1..5] do for j in [i+1..5] do
  M3lift[i][j] := NaiveLiftPol(M3[i][j], K);
  M3lift[j][i] := -M3lift[i][j];
end for; end for;
M4lift := ZeroMatrix(Rlift, 5, 5);
for i in [1..5] do for j in [i+1..5] do
  M4lift[i][j] := NaiveLiftPol(M4[i][j], K);
  M4lift[j][i] := -M4lift[i][j];
end for; end for;

Rlift<x> := PolynomialRing(K);
Slift<y> := PolynomialRing(Rlift);
SSlift<y1,y2,y3> := PolynomialRing(Rlift,3);
philift := hom<SSlift -> Slift | 0, 0, y>;

Mlift := ZeroMatrix(SSlift,5,5);
Mlift := M1lift + M2lift*y3 + M3lift*y2 + M4lift*y1;

print "Lifted to characteristic zero, now projecting to plane";


flift := philift(Basis(EliminationIdeal(Ideal(Pfaffians(Mlift, 4)),2))[1]);

print "Clearing denominators...";
// clear denominators;
listofdenoms := [];
for c in Coefficients(flift) do
  if Degree(F) eq 1 then
    listofdenoms cat:=[Denominator(cc) : cc in Coefficients(c)];
  else
    for cc in Coefficients(c) do
      listofdenoms cat:=[Denominator(e) : e in Eltseq(cc)];
    end for;
  end if;
end for;
flift *:= LCM(listofdenoms);

print "Making monic...";

flift := y^5 + Coefficient(flift,4)*y^4
         + Coefficient(flift,3)*Coefficient(flift,5)*y^3
         + Coefficient(flift,2)*Coefficient(flift,5)^2*y^2
         + Coefficient(flift,1)*Coefficient(flift,5)^3*y
         + Coefficient(flift,0)*Coefficient(flift,5)^4; // MAKE MONIC

print "----";

// PRINT TO FILE "inputcurve_4.m", CAN BE FED TO TUITMAN'S ALGORITHM
// (IF F IS A PRIME FIELD, FEED TO pcc_p.m, ELSE FEED TO pcc_q.m)

if Degree(F) eq 1 then
  Write("inputcurve_5.m", "p := "*IntegerToString(#F)*"; \n" : Overwrite := true);
else
  Write("inputcurve_5.m", "q := "*IntegerToString(#F)*"; \n" : Overwrite := true);
end if;

Write("inputcurve_5.m", "Q := \n");
Write("inputcurve_5.m", flift);
Write("inputcurve_5.m", ";");
