// This Magma script should compute a non-singular model of the
// genus 5 and genus 7 curves. It involves the potentially length computation
// of the canonical model, but I don't see a way around that.

load "gl2data.txt";

// The function below uses the models and maps we have computed to
// map points down to the j-line. It takes two parameters -
// a subgroup number, and a list of rationals.

function maptoj(subnum,ptlist);
  curgp := subnum;
  curptlist := ptlist;
  while (curgp gt 1) do
    covergp := newsublist[curgp][7];
    mp := eval Read(Sprintf("%omap%o.txt",curgp,covergp));  
    curpt := Domain(mp)!curptlist;
    newpt := mp(curpt);
    curptlist := Coordinates(newpt);
    curgp := covergp;
  end while;
  return newpt;
end function;

subnumlist := [674,675,676,677,678,679,680,683,687,688,690,691,692,693,694,695,706,712,713,714];
Q := Rationals();
P1Q := ProjectiveSpace(Q,1);
R<w,x,y,z> := PolynomialRing(Q,4);
RR<t> := PolynomialRing(Q);
jlist := [ 
[P1Q![54000,1]],
[P1Q![1,0],P1Q![16581375,1]],
[P1Q![0,1]],
[P1Q![1,0],P1Q![16581375,1]],
[P1Q![1,0],P1Q![16581375,1]],
[P1Q![0,1]],
[P1Q![54000,1]],
[P1Q![54000,1]],
[P1Q![0,1]],
[P1Q![0,1]],
[P1Q![1,0],P1Q![16581375,1]],
[P1Q![1,0],P1Q![-3375,1]],
[P1Q![1,0],P1Q![-3375,1]],
[P1Q![1,0],P1Q![-3375,1]],
[P1Q![1,0],P1Q![-3375,1]],
[P1Q![54000,1]],
[P1Q![1,0],P1Q![1728,1]],
[P1Q![1,0],P1Q![-3375,1]],
[P1Q![1,0],P1Q![16581375,1]],
[P1Q![0,1]]
];

letlist := ["a","b","c","d","e","f","g"];
for ind in [1..#subnumlist] do
  subnum := subnumlist[ind];
  covnum := newsublist[subnum][7];

  // Compute and parse double cover of elliptic curve

  C := eval Read("eq" cat IntegerToString(subnum) cat ".txt");
  E := eval Read("eq" cat IntegerToString(covnum) cat ".txt");
  gen := Genus(C);

  de := DefiningEquations(C);
  polylist := [ Evaluate(de[i],[x,y,1,w]) : i in [1..#de]];
  I := ideal<R | polylist>;
  G := GroebnerBasis(I);
  poly := G[1];
  cofs := Coefficients(poly);
  denom := LCM([ Denominator(cofs[i]) : i in [1..#cofs]]);
  poly := poly*denom;
  doublecover := -Evaluate(poly,[0,x,y,1]);

  // Step 1 - Compute canonical image

  printf "Computing canonical map for X%o.\n",subnum;
  phi := CanonicalMap(C);
  printf "Done. Computing canonical image.\n";
  D := CanonicalImage(Domain(phi),phi);
  AssignNames(~D,[ letlist[i] : i in [1..gen]]);
  printf "The canonical image is %o.\n",D;
  printf "Done. Restricting map.\n";
  phi := Restriction(phi,C,D);

  // Step 2 - Write down points on C with x,y,z in Q and w in a quadratic field
  // Compute their images on D, and their images under the automorphism on D.

  printf "Creating points on C over quadratic fields and mapping them to the canonical model.\n";
  G, MWmap := MordellWeilGroup(E);
  ptlist := <>;
  done := false;
  cur := 0;
  while done eq false do
    cur := cur + 1;
    ptonE := MWmap(cur*G.2);
    xval := Q!ptonE[1];
    yval := Q!ptonE[2];
    zval := Q!1;
    wvalsqr := Evaluate(doublecover,[0,xval,yval,0]);
    if not IsSquare(wvalsqr) then
      K<ww> := NumberField(t^2 - wvalsqr);
      ptonC := C(K)![xval,yval,1,ww];
      ptonC2 := C(K)![xval,yval,1,-ww];
  // Step 3 - compute their images on D
      Append(~ptlist,<ptonE,phi(ptonC),phi(ptonC2)>);
      printf "Created %o point(s) so far.\n",#ptlist;
    end if;
    if #ptlist ge 2*gen then
      done := true;
    end if;
  end while;

  // Step 4 - Use these points to guess the automorphism on D

  M := ZeroMatrix(Rationals(),0,gen^2);
  done := false;
  it := 0;
  while done eq false do
    it := it + 1;
    pt := ptlist[it];
    for j in [2..gen] do
      vec := ZeroMatrix(Rationals(),2,gen^2);
      for k in [1..gen] do
        e1 := pt[2][k]*pt[3][j];
        e2 := pt[2][k]*pt[3][1];
        if Degree(MinimalPolynomial(e1)) eq 1 then
          E1 := [ e1, 0 ];
        else
          E1 := Eltseq(e1);
        end if;
        if Degree(MinimalPolynomial(e2)) eq 1 then
          E2 := [ e2, 0 ];
        else
          E2 := Eltseq(e2);
        end if;
        vec[1][k] := E1[1];
        vec[1][(j-1)*gen+k] := -E2[1];
        vec[2][k] := E1[2];
        vec[2][(j-1)*gen+k] := -E2[2];
      end for;
      M := VerticalJoin(M,vec);
    end for;
    if Rank(M) eq (gen^2 - 1) then
      done := true;
    end if;
  end while;
  N := Nullspace(Transpose(M));
  v := N.1;
  elt := Eltseq(v);
  m := Matrix(Rationals(),gen,gen,elt);
  printf "Matrix of automorphism is %o.\n",m;
  // Check that this process worked
  vars := [ D.i : i in [1..gen]];
  maplist := [ &+[ m[i][j]*vars[j] : j in [1..gen]] : i in [1..gen]];
  aut := iso<D->D|maplist,maplist>;

  // Step 5 - create a map to P^(g-2) by omitting the variable
  // that is inverted by the automorphism.
  if IsDiagonal(m) then
    if #[i : i in [1..gen] | m[i][i] eq -1] eq 1 then
      badvar := Index([ m[i][i] : i in [1..gen]],-1);
      lst := [ D.j : j in [1..gen] | j ne badvar];
      quomap := map<D->ProjectiveSpace(Rationals(),gen-2) | lst>;
      X := Image(quomap);
      quomap := Restriction(quomap,D,X);
      printf "The curve X is %o.\n",X;
      ptonX := X(Q)!quomap(ptlist[1][2]);
      E1, emap := EllipticCurve(X,ptonX);
      emap := Extend(emap);
      a, emapinv := IsInvertible(emap);
      //emapinv := Extend(emapinv);
      printf "The curve E1 is %o.\n",E1;
      E2, emap2 := MinimalModel(E1);
      printf "The curve E2 is %o.\n",E2;
      if not (aInvariants(E) eq aInvariants(E2)) then
        printf "Error. We didn't get exactly the same minimal model.\n";
        bad := 0;
        bad2 := 1/bad;
      end if;
      bigmap := quomap*emap*emap2;
    else
      printf "The automorphism doesn't have the right shape!\n";
      bad := 0;
      bad2 := 1/bad;
    end if;
  else
    printf "The automorphism doesn't have the right shape!\n";
    bad := 0;
    bad2 := 1/bad;
  end if;

  // Step 6 - Elliptic curves with positive rank have infinitely many
  // automorphisms. Find the right automorphism that gives the exact map
  // from D -> E2. We can take quomap*emap*emap2 and then post-compose
  // with a translation and/or [-1] to get things to match.
  P1 := E2!ptlist[2][1];
  P2 := E2!ptlist[3][1];
  Q1 := E2(Q)!bigmap(ptlist[2][2]);
  Q2 := E2(Q)!bigmap(ptlist[3][2]);
  if P2 eq (Q2 + P1 - Q1) then
    // The automorphism is psi(R) = R + (P1-Q1)
    emap3 := TranslationMap(E2,P1-Q1);
  end if;
  if P2 eq (-Q2 + P1 + Q1) then
    // The automorphism is psi(R) = -R + (Q1 + P1)
    emap3 := NegationMap(E2)*TranslationMap(E2,Q1+P1);
  end if;
  // This is a map, but not a map of schemes. That's OK though.
  printf "The map from the canonical model of X%o to X%o is the following.\n",subnum,covnum;
  bigmap2 := bigmap*emap3;
  printf "%o.\n",bigmap2;

  // Step 7 - Compute preimage on Xcov of the given points on the jline
  done := false;
  gplist := [covnum];
  while done eq false do
    cur := gplist[#gplist];
    cover := newsublist[cur][7];
    if (cover gt 0) then
      Append(~gplist,cover);
    else 
      done := true;
    end if;
  end while;
  Reverse(~gplist);

  ptlist := jlist[ind];
  printf "Computing preimages of j-invariants on X%o.\n",covnum;
  for i in [1..#gplist-1] do
    covmap := eval Read(IntegerToString(gplist[i+1]) cat "map" cat IntegerToString(gplist[i]) cat ".txt");
    covmap := Extend(covmap);
    curspace := Codomain(covmap);
    covspace := Domain(covmap);
    newptlist := [ ];
    for j in [1..#ptlist] do
      X := (curspace!Eltseq(ptlist[j]))@@covmap;
      lst := RationalPoints(X);
      for p in lst do
        Append(~newptlist,covspace!p);
      end for;
    end for;
    ptlist := newptlist;
  end for;  
  printf "Done.\n";

  // Step 8 - Compute all the rational points on the genus 5 or 7 curve.    
  ptlist2 := [];
  bigmap2 := quomap*emap*emap2*emap3;
  for p in ptlist do
    q := p@@emap3;
    r := q@@emap2;
    s := emapinv(r);
    X := s@@quomap;
    S := RationalPoints(X);
    for s in S do
      Append(~ptlist2,<D!Eltseq(s),P1Q!Eltseq(maptoj(covnum,Eltseq(bigmap2(Domain(bigmap2)!Eltseq(s)))))>);
    end for;    
  end for;
  ptlist3 := [];
  for i in [1..#ptlist2] do
    found := false;
    for j in [1..i-1] do
      if ptlist2[i] eq ptlist2[j] then
        found := true;
      end if;
    end for;
    if (found eq false) then
      Append(~ptlist3,ptlist2[i]);
    end if;
  end for;
  printf "There are %o points on X%o.\n",#ptlist3,subnum;
  printf "The only points on X%o are the following: %o.\n",subnum,ptlist3;
end for;