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

// PARAMETERS

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

print "----";
print "Testing our lifting procedure for a generic tetragonal curve of genus", g;
print "over GF(", Characteristic(F),",", Degree(F), "), with Maroni invariants", e[1], e[2], e[3];
print "and Schreyer invariants", b[1], b[2];
print "The Maroni invariants should sum up to g - 3, which is", &+e eq g-3;
print "The Schreyer invariants should sum up to g - 5, which is", &+b eq g-5;
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 4 COVER WITH ABOVE PARAMETERS 
// AND SCRAMBLE OUTPUT SOMEWHAT

RandomQuadric := function(F,e,b)
  R<x> := PolynomialRing(F);
  S<y1,y2> := PolynomialRing(R,2);
  quad := S ! 0;
  if #[0..e[3] + e[3] - b] ge 1 then quad +:= &+[Random(F)*x^i : i in [0..e[3] + e[3] - b]]; end if;
  if #[0..e[1] + e[3] - b] ge 1 then quad +:= &+[Random(F)*x^i : i in [0..e[1] + e[3] - b]]*y1; end if;
  if #[0..e[2] + e[3] - b] ge 1 then quad +:= &+[Random(F)*x^i : i in [0..e[2] + e[3] - b]]*y2; end if;
  if #[0..e[1] + e[1] - b] ge 1 then quad +:= &+[Random(F)*x^i : i in [0..e[1] + e[1] - b]]*y1^2; end if;
  if #[0..e[1] + e[2] - b] ge 1 then quad +:= &+[Random(F)*x^i : i in [0..e[1] + e[2] - b]]*y1*y2; end if;
  if #[0..e[2] + e[2] - b] ge 1 then quad +:= &+[Random(F)*x^i : i in [0..e[2] + e[2] - b]]*y2^2; end if;
  return quad;
end function;

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

Q1 := SS ! RandomQuadric(F, e, b[1]);
Q2 := SS ! RandomQuadric(F, e, b[2]);
f := phi(Resultant(Q1, Q2, 2));
f := y^4 + Coefficient(f,3)*y^3 
     + Coefficient(f,2)*Coefficient(f,4)*y^2
     + Coefficient(f,1)*Coefficient(f,4)^2*y
     + Coefficient(f,0)*Coefficient(f,4)^3; // 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 TWO QUADRATIC FORMS

BhargavaQuadrics := function(f)
  kx := BaseRing(Parent(f)); frkx := FieldOfFractions(kx);
  Mkx := MatrixAlgebra(kx, 3);
  f0 := Coefficients(f)[1]; f1 := Coefficients(f)[2]; f2 := Coefficients(f)[3]; f3 := Coefficients(f)[4];
  fres := Parent(f) !  [-f3^2*f0 + 4*f0*f2 - f1^2, f1*f3 - 4*f0, -f2, 1];
  FF<Y> := FunctionField(f); FFres<Yres> := FunctionField(fres);
  print "Genus of resolvent curve should be g + 1 which is", Genus(FFres) eq gf + 1;

  // COMPUTE REDUCED BASIS AND BASE CHANGES TO NAIVE BASES
  red_bas := ShortBasis(DivisorGroup(FF) ! 0);
  n_bas := [1, -f0*Y^(-1), Y, Y^2];
  red_bas_res := ShortBasis(DivisorGroup(FFres) !  0);
  n_bas_res := [1, Yres, Yres^2 - f2*Yres + f1*f3-4*f0];

  N, h := Module(n_bas, frkx :  IsBasis := true);
  Mres, hres := Module(red_bas_res, kx :  IsBasis := true);
  A := Matrix(frkx, 3, 3,
              [Inverse(h)(red_bas[i])[j] :  i in [2,3,4], j in [2,3,4]]);
  B := Matrix(frkx, 2, 2,
              [Inverse(hres)(n_bas_res[i])[j] :  i in [2,3], j in [2,3]]);

  // COMPUTE TERNARY QUADRATIC FORMS
  Q1 := Matrix(kx, 3, 3, [f0, 0, f1/2, 0, 1, -f3/2, f1/2, -f3/2 ,f2]);
  Q2 := Matrix(kx, 3, 3, [0, kx ! -1/2, f3/2, kx ! -1/2, 0, 0, f3/2, 0, 1]);
  Q1 := Transpose(A)*Q1*A; Q2 := Transpose(A)*Q2*A;
  return Mkx ! (B[1,1]*Q1 + B[1,2]*Q2), Mkx ! (B[2,1]*Q1 + B[2,2]*Q2);
end function;

print "Computing Bhargava's quadrics...";
//tim := Cputime();
Q1, Q2 := BhargavaQuadrics(f);
//print "Took time", Cputime(tim);

Rlift<x> := PolynomialRing(K);
Q1lift := Matrix(Rlift, 3, 3, 
             [NaiveLiftPol(Q1[j][i], K) : i in [1,2,3], j in [1,2,3]]);
Q2lift := Matrix(Rlift, 3, 3,
             [NaiveLiftPol(Q2[j][i], K) : i in [1,2,3], j in [1,2,3]]);

// PROJECT DOWN AND MAKE MONIC

MatrixToQuadric := function(Q)
  R := BaseRing(Q);
  S<y1,y2> := PolynomialRing(R,2);
  return Q[1][1]*y1^2 + 2*Q[1][2]*y1*y2 + 2*Q[1][3]*y1 + Q[2][2]*y2^2 + 2*Q[2][3]*y2 + Q[3][3];
end function;

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

f1lift := SSlift ! MatrixToQuadric(Q1lift);
f2lift := SSlift ! MatrixToQuadric(Q2lift);

print "Computing resultant...";
//tim := Cputime();
flift := philift(Resultant(f1lift,f2lift,1));
//print "Took time", Cputime(tim);

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

// 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)

print "----";

if Degree(F) eq 1 then
  Write("inputcurve_4.m", "p := "*IntegerToString(#F)*"; \n" : Overwrite := true);
  print "Writing output to inputcurve_4.m, can be fed as input to pcc_p";
else
  Write("inputcurve_4.m", "q := "*IntegerToString(#F)*"; \n" : Overwrite := true);
  print "Writing output to inputcurve_4.m, can be fed as input to pcc_q";
end if;

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