// This Magma script is going to be used to combine all the data we have
// to test an input elliptic curve, and output what the image of the
// 2-adic Galois representation is.

load "gl2data.txt";
load "gl2finedata.txt";

F<t> := FunctionField(Rationals());
R<x> := PolynomialRing(Rationals());

function fileexists(str)
  retval := System("test -e " cat str);
  retval := retval div 256;
  ret := false;
  if (retval eq 0) then
    ret := true;
  end if;
  return ret;
end function;

// Given a genus zero curve X_G with chosen function field generator 
// Q(X_G) = Q(t), return j as a rational function in t.

function jmap(subnum)
  curmap := F!t;
  done := false;
  curgp := subnum;
  while (done eq false) do
    if (curgp eq 1) then
      done := true;
      j := curmap;
    else
      covernum := newsublist[curgp][7];
      filename := Sprintf("%omap%o.txt",curgp,covernum);
      maptocover := eval Read(filename);
      newcurmap := maptocover([ curmap, 1])[1];
      curmap := newcurmap;
      curgp := covernum;
    end if;
  end while;
  return curmap;
end function; 

// Hard-coded list of j-invariants for genus 2, 3, 5 and 7 subgroups.

genustwoorhigher := [
<441,-18234932071051198464000/48661191875666868481>,
<441,-35817550197738955933474532061609984000/2301619141096101839813550846721>,
<556,78608>,
<558,68769820673/16>,
<563,2048>,
<566,16974593/256>,
<619,-631595585199146625/218340105584896>,
<649,777228872334890625/60523872256>];

// Rank one elliptic curves. All the genus one curves either have
// no local points, rank zero and all torsion points are cusps/CM, or
// rank 1. All covers of rank one elliptic curves have genus > 1.

rankoneell := [150, 153, 155, 156, 165, 166, 167, 281, 284, 288, 289, 291, 295, 297, 302, 304, 308, 309, 312, 318, 320, 323, 324, 326, 328, 349, 350];

// Genus zero curves that have points

genzeropoints := [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,15, 16, 17, 
18, 19, 20, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 
39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 55, 56, 57, 58, 60, 61, 62, 63, 
64, 65, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 
85, 86, 87, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104,
105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 
121, 122, 123, 124, 177, 178, 181, 183, 185, 187, 188, 189, 190, 191, 192, 193, 
194, 195, 196, 197, 199, 200, 202, 203, 204, 205, 207, 208, 209, 210, 211, 212, 
213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 
229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243 ];

printf "Enter the a-invariants of an elliptic curve:\n";
printf "Enter xy coeff:";
readi a1;
printf "Enter x^2 coeff:";
readi a2;
printf "Enter y coeff:";
readi a3;
printf "Enter x coeff:";
readi a4;
printf "Enter constant coeff:";
readi a6;
ainvars := [a1,a2,a3,a4,a6];
E := EllipticCurve(ainvars);
E := MinimalModel(E);
printf "Curve has conductor %o.\n",Conductor(E);

done1 := false;
curgp := 1;
// Make a list of subgroups, and all the choices for the invariants
invarlist := [];
jval := jInvariant(E);
printf "j-invariant is %o.\n",jval;
P1 := ProjectiveSpace(Rationals(),1);
Append(~invarlist,<1,[P1![jval,1]]>);

while done1 eq false do
  // For each maximal subgroup G of curgp, test to see if the points on X_curgp
  // lift to points on X_G.
  testmore := false;
  moregps := [];
  printf "Current group is %o.\n",curgp;
  for ind in newsublistmaxsub[curgp] do
    printf "Working on X_%o.\n",ind;
    if newsublist[ind][5] eq 0 then
      if ind in genzeropoints then
        jfunc := jmap(ind);
        num := R!Numerator(jfunc);
        denom := R!Denominator(jfunc);
        poly := num - jval*denom;
        rtlist := Roots(poly);
        p1points := [ P1![ rtlist[i][1],1] : i in [1..#rtlist]];
        if Degree(num) eq Degree(denom) then
          if LeadingCoefficient(num)/LeadingCoefficient(denom) eq jval then
            Append(~p1points, P1![1,0]);
          end if;
        end if;
        if #p1points ge 1 then
          printf "Point on X_%o lifts to a point on X_%o!\n",curgp,ind;
          testmore := true;
          Append(~invarlist,<ind,p1points>);
          Append(~moregps,ind);
        end if;                   
      end if;
    end if;
    if newsublist[ind][5] eq 1 then
      if ind in rankoneell then
        // Is the cover of this group in the invarlist?
        ind2 := Index([ invarlist[i][1] : i in [1..#invarlist]],newsublist[ind][7]);
        if (ind2 ne 0) then
          p1points := invarlist[ind2][2];
        else
          coverind := newsublist[ind][7];
          jfunc := jmap(coverind);
          num := R!Numerator(jfunc);
          denom := R!Denominator(jfunc);
          poly := num - jval*denom;
          rtlist := Roots(poly);
          p1points := [ P1![ rtlist[i][1],1] : i in [1..#rtlist]];
          if Degree(num) eq Degree(denom) then
            if LeadingCoefficient(num)/LeadingCoefficient(denom) eq jval then
              Append(~p1points, P1![1,0]);
            end if;
          end if;
        end if;
        ellmap := eval Read(IntegerToString(ind) cat "map" cat IntegerToString(newsublist[ind][7]) cat ".txt");
        ptset := { Codomain(ellmap)!Eltseq(p1points[i]) : i in [1..#p1points]};
        //printf "Computing preimage of map from X%o to X%o.\n",ind,newsublist[ind][7];
        //printf "Point set is %o.\n",ptset;
        addinvarlist := <>;
        for p in ptset do
          preims := RationalPoints(p@@ellmap);
          //printf "The preimages of %o are %o.\n",p,preims;
          for p2 in preims do
            //printf "The image of %o is %o.\n",p2,ellmap(p2);
            if ellmap(p2) eq p then
              Append(~addinvarlist,p2);
            end if;
          end for;
        end for;
        if #addinvarlist gt 0 then
          printf "Point on X_%o lifts to a point on X_%o!\n",curgp,ind;
          testmore := true;
          Append(~moregps,ind);
        end if;                 
      end if;       
    end if;
    if newsublist[ind][5] ge 2 then
      ind2 := Index( [ genustwoorhigher[i][2] : i in [1..#genustwoorhigher]],jval);
      if (ind2 ne 0) then
        if ind eq genustwoorhigher[ind2][1] then
          printf "Point on X_%o lifts to a point on X_%o!\n",curgp,ind;
          testmore := true; 
          Append(~moregps,ind);       
        end if;
      end if;
    end if;
  end for;
  if #moregps gt 0 then
    curgp := moregps[#moregps];
  end if;
  if testmore eq false then
    done1 := true;
    minusonegp := curgp;
  end if;
end while;
printf "The image for the subgroup containing -I is X%o.\n",minusonegp;
// Deal with subgroups in finesublist
finesubs := [ i : i in [1..#finesublist] | finesublist[i][1] eq minusonegp];
index := newsublist[minusonegp][1];
subname := "X" cat IntegerToString(minusonegp);
image := newsublist[minusonegp][3];
foundtimes := 0;
for n in [1..#finesubs] do
  foundhere := false;
  // If we get here, the genus of X_minusonegp is zero.
  ind := Index([ invarlist[i][1] : i in [1..#invarlist]],minusonegp);
  tvals1 := invarlist[ind][2];
  tvals2 := [];
  infinitetval := [];
  for m in [1..#tvals1] do
    if tvals1[m][2] eq 0 then
      // Point at infinity.
      Append(~infinitetval,true);
      Append(~tvals2,0);
    else
      Append(~tvals2,Rationals()!tvals1[m][1]);
      Append(~infinitetval,false);
    end if;
  end for;
  fileopenname := "eqf" cat IntegerToString(finesubs[n]) cat ".txt";
  // Check to see if file exists yet
  if fileexists(fileopenname) then
    Et := eval Read("eqf" cat IntegerToString(finesubs[n]) cat ".txt");
    ainvars := aInvariants(Et);
    A := ainvars[4];
    B := ainvars[5];
    for m in [1..#tvals2] do
      if infinitetval[m] eq false then
        Aspec := Evaluate(A,tvals2[m]);
        Bspec := Evaluate(B,tvals2[m]);
        E2 := EllipticCurve([0,0,0,Aspec,Bspec]);
        chk := IsIsomorphic(E,E2);
        if (chk eq true) then
          foundhere := true;
        end if;
      else
        //printf "Extending cover to Et -> P^1.\n";
        //printf "A = %o.\n",A;
        //printf "B = %o.\n",B;
        degA := Degree(Numerator(A));
        degB := Degree(Numerator(B));
        powt := Max(Ceiling(degA/4),Ceiling(degB/6));
        A2 := Evaluate(A,1/t)*t^(4*powt);
        B2 := Evaluate(B,1/t)*t^(6*powt);
        Aspec := Evaluate(A2,0);
        Bspec := Evaluate(B2,0);
        //printf "Curve is y^2 = x^3 + %ox + %o.\n",Aspec,Bspec;
        if (-64*Aspec^3 - 432*Bspec^2) ne 0 then
          E2 := EllipticCurve([0,0,0,Aspec,Bspec]);
          chk := IsIsomorphic(E,E2);
          if (chk eq true) then
            foundhere := true;
          end if;
        end if;
      end if;
    end for;
  else
    printf "Sorry, we can't check %o. The model is not available.\n",finesublist[finesubs[n]][5];
  end if;
  if foundhere eq true then
    foundtimes := foundtimes + 1;
    subname := finesublist[finesubs[n]][5];
    //printf "The image should be %o.\n",subname;
    image := finesublist[finesubs[n]][4];
    index := finesublist[finesubs[n]][2];
  end if;
end for;
if foundtimes gt 1 then
  printf "Error! The elliptic curve was found multiple times!\n";
  printf "Something is haywire!\n";
  bad := 0;
  bad2 := 1/bad;
end if;
printf "Index of 2-adic image: %o.\n",index;
printf "Name of curve that E lies on: %o.\n",subname;