// This Magma script is being used to compute the resolvent polynomials
// that the Dokchitser's wrote down for use in determining the
// image of elliptic curve Galois representations. Right now,
// it tries to prove that the Galois representation attached to a particular
// elliptic curve surjects onto the group it is supposed to.

load "gl2data.txt";

function liftsub(G,n,m)
  H2 := GL(2,Integers(2^m));
  H3 := GL(2,Integers(2^n));
  genlist := [];
  for g in Generators(G) do
    Append(~genlist,H2!g);
  end for;
  phi := hom<H2->H3 | [H3!H2.i : i in [1..#Generators(H2)]]>;
  ker := Kernel(phi);
  // Add generators for the kernel of GL_2(Z/2^m Z) -> GL_2(Z/2^n Z)
  return sub<H2 | genlist, ker>;
end function;

// E := EllipticCurve([0,0,1,-3014658660,150916472601529]);
// E := EllipticCurve([0,0,1,-26668521183591560,-1676372876256754456631251]);
// subnum := 441;
E := EllipticCurve([0,1,0,-28,48]);
subnum := 556;
// E := EllipticCurve([1,0,1,-3500,-79974]);
// subnum := 558;
// E := EllipticCurve([0,1,0,-3,-2]);
// subnum := 563;
// E := EllipticCurve([1,0,1,-220,-1254]);
// subnum := 566;
// E := EllipticCurve([1,-1,1,-303870,-81409651]);
// subnum := 619;
// E := EllipticCurve([1,-1,1,-325630,-71434867]);
// subnum := 649;

printf "Working on curve %o.\n",aInvariants(E);
printf "Testing subgroup %o.\n",subnum;

// Assume that m is not equal to 2.

rightval := Max([ newsublist[subnum][4][i][2] : i in [1..#newsublist[subnum][4]]]);
G := newsublist[subnum][3];
if #BaseRing(G) lt 2^(rightval) then
  printf "Increasing modulus from %o to %o.\n",#BaseRing(G),2^rightval;
  G := liftsub(G,Valuation(#BaseRing(G),2),rightval);
end if;
m := 2^rightval;
conjset := Conjugates(GL(2,Integers(m)),G);
seq := SetToSequence(conjset);

// Which is the right conjugate of G?

// Use the elliptic function x+y
// Use h(x) = x^2

// Go over all points of order m
eltorderm := [ [i,j] : i in [1..m], j in [1..m] | GCD([i,j,m]) eq 1 ];

T<x> := PolynomialRing(Integers());

if (m eq 4) then
  decprec := 100;
  decprec2 := 100;
  h := x^2;
end if;
if (m eq 8) then
  decprec := 1000;
  decprec2 := 1000;
  h := x^3;
end if;
if (m eq 16) then
  decprec := 1000;
  decprec2 := 1000;
  h := x^3;
end if;
if (subnum eq 441) then
  decprec := 10000;
  decprec2 := 6000;
  h := x^3;
end if;
if (subnum eq 556) then
  decprec := 1000;
  decprec2 := 500;
  h := x^3;
end if;
if (subnum eq 558) then
  decprec := 1000;
  decprec2 := 500;
  h := x^3;
end if;
if (subnum eq 563) then
  decprec := 1000;
  decprec2 := 500;
  h := x^3;
end if;
if (subnum eq 566) then
  decprec := 1000;
  decprec2 := 500;
  h := x^3;
end if;
if (subnum eq 619) then
  decprec := 4000;
  decprec2 := 1000;
  h := x^3;
end if;
if (subnum eq 649) then
  decprec := 4000;
  decprec2 := 2000;
  h := x^3;
end if;

polylist := [];
C<I> := ComplexField(decprec);
CCC<I2> := ComplexField(decprec2);
R<z> := PolynomialRing(C);
periodlist := Periods(E : Precision := decprec);
Rsmall := RealField(50);
RSTUV<zsmall> := PolynomialRing(Rsmall);

printf "Enumerating points of order %o.\n",m;
hvals1 := [];
ellvals1 := [];
hvals2 := [];
ellvals2 := [];
count := 0;
for pt in eltorderm do
  count := count + 1;
  if (count mod 20 eq 0) then
    printf "Done %o of %o.\n",count,#eltorderm;
  end if;
  ellpt := EllipticExponential(E,periodlist[1]*pt[1]/m+periodlist[2]*pt[2]/m);
  ellval := 2*ellpt[1] + 4*ellpt[2];
  // ellval is an algebraic integer according to Theorem VII.3.4 from Silverman
  Append(~hvals1,Evaluate(h,ellval));
  Append(~ellvals1,ellval);
  Append(~hvals2,CCC!Evaluate(h,ellval));
  Append(~ellvals2,CCC!ellval);
end for;

// Figure out which conjugate we should pick.

printf "Selecting which conjugate of G we should use.\n";
best := 0;
minerr := -1;
for choice in [1..#seq] do
  printf "Testing choice %o. ",choice;
  candidate := seq[choice];
  num := C!0;
  pt1 := [1,m];
  pt2 := [m,1];
  pt3 := [1,1];
  for g in candidate do
    a := pt1[1];
    b := pt1[2];
    impt := [ g[1][1]*a + g[2][1]*b, g[1][2]*a + g[2][2]*b ];
    impt := [ Integers()!impt[1], Integers()!impt[2]];
    if impt[1] eq 0 then
      impt[1] := m;
    end if;
    if impt[2] eq 0 then
      impt[2] := m;
    end if;
    ind1 := Index(eltorderm,pt1);
    ind2 := Index(eltorderm,impt);
    term1 := hvals2[ind2];
    a := pt2[1];
    b := pt2[2];
    impt := [ g[1][1]*a + g[2][1]*b, g[1][2]*a + g[2][2]*b ];
    impt := [ Integers()!impt[1], Integers()!impt[2]];
    if impt[1] eq 0 then
      impt[1] := m;
    end if;
    if impt[2] eq 0 then
      impt[2] := m;
    end if;
    ind1 := Index(eltorderm,pt2);
    ind2 := Index(eltorderm,impt);
    term2 := hvals2[ind2];
    a := pt3[1];
    b := pt3[2];
    impt := [ g[1][1]*a + g[2][1]*b, g[1][2]*a + g[2][2]*b ];
    impt := [ Integers()!impt[1], Integers()!impt[2]];
    if impt[1] eq 0 then
      impt[1] := m;
    end if;
    if impt[2] eq 0 then
      impt[2] := m;
    end if;
    ind1 := Index(eltorderm,pt3);
    ind2 := Index(eltorderm,impt);
    term3 := hvals2[ind2];
    num := num + (term1*term2 + term1*term3 + term2*term3);
  end for;
  err := RealField(10)!AbsoluteValue(num - Round(Real(num)));
  printf "We get error %o.\n",err;
  if (minerr eq -1) or (err lt minerr) then
    best := choice;
    minerr := err;
  end if;  
end for; 
printf "Choice %o is best.\n",best;
G := seq[best];
printf "The order of G is %o.\n",#G;
CC := ConjugacyClasses(G);
printf "There are %o conjugacy classes of G.\n",#CC;
printf "Their sizes are %o.\n",[ CC[i][2] : i in [1..#CC]];

printf "Computing fpolynomial.\n";
fpolycomp := &*[ z - ellvals1[i] : i in [1..#ellvals1]];
fpolycheck := &*[ z + (Rsmall!AbsoluteValue(ellvals1[i])) : i in [1..#ellvals1]]
;
newcoefflist := [ Round(Real(Coefficient(fpolycomp,i))) : i in [0..Degree(fpolycomp)]];
maxerr := Max([ RealField(10)!AbsoluteValue(Coefficient(fpolycomp,i) - newcoefflist[i+1]) : i in [0..Degree(fpolycomp)]]);
fpoly := &+[ newcoefflist[i]*x^(i-1) : i in [1..#newcoefflist]];
printf "Maximal error on fpolynomial was %o.\n",maxerr;
maxcofsize := Max([ AbsoluteValue(Coefficient(fpolycheck,i)) : i in 
[0..Degree(fpolycheck)]]);
printf "Maximal coefficient size of f polynomial was %o.\n",
RealField(10)!maxcofsize;
if (maxcofsize gt 10^decprec) then
  printf "The coefficients are too big. We need more decimal precision!\n";
  bad := 0;
  bad2 := 1/bad;
end if;

classlist := [];
for i in [1..#CC] do
  // Create polynomial for conjugacy class i
  printf "Working on class %o. ",i;
  cent := Centralizer(G,CC[i][3]);
  class := [ s^(-1)*CC[i][3]*s : s in Transversal(G,cent) ];
  Append(~classlist,class);
  poly := R!1;
  polycheck := RSTUV!1;
  for g in class do
    term := C!0;
    for pt in eltorderm do
      a := pt[1];
      b := pt[2];
      // The action is on the right, just like models9.txt
      impt := [ g[1][1]*a + g[2][1]*b, g[1][2]*a + g[2][2]*b ];
      impt := [ Integers()!impt[1], Integers()!impt[2]];
      if impt[1] eq 0 then
        impt[1] := m;
      end if;
      if impt[2] eq 0 then
        impt[2] := m;
      end if;
      ind1 := Index(eltorderm,pt);
      ind2 := Index(eltorderm,impt);
      term := term + hvals2[ind1]*ellvals2[ind2];      
    end for;
    poly := poly*(z-term);
    polycheck := polycheck*(z + (Rsmall!AbsoluteValue(term)));
  end for;
  roundpoly := 0;
  maxcurerror := RealField(10)!0;
  for n in [0..Degree(poly)] do
    newcoeff := Round(Real(Coefficient(poly,n)));
    rounderr := RealField(10)!AbsoluteValue(Coefficient(poly,n)-newcoeff);
    maxcurerror := Max(rounderr,maxcurerror);
    roundpoly := roundpoly + newcoeff*x^n;
  end for;
  maxcoeff := Max([ AbsoluteValue(Coefficient(polycheck,ijk)) : ijk in 
  [0..Degree(polycheck)]]);
  maxerr := Max(maxerr,maxcurerror);  
  Append(~polylist,roundpoly);    
  printf "Maximal error on polynomial was %o.\n",maxcurerror;
  if (maxcoeff gt 10^decprec2) then
    printf "The coefficient size is too big. We need more decimal precision.\n";
    bad := 0;
    bad2 := 1/bad;
  end if;
end for;
printf "Maximal error was %o.\n",maxerr;

// GCD check
for i in [1..#polylist] do
  for j in [i+1..#polylist] do
    deg := Degree(GCD(polylist[i],polylist[j]));
    if (deg gt 0) then
      printf "Error! Polynomials %o and %o have a common factor.\n",i,j;
    end if;
  end for;
end for;

printf "Starting the Frobenius computations.\n";
classesfound := {};
primecount := 0;
for p in [2..30000] do
  if IsPrime(p) then
    primecount := primecount + 1;
    if (primecount mod 1000 eq 0) then
      printf "Tested %o primes.\n",primecount;
    end if;
    if GCD(p,2*Conductor(E)) eq 1 then
      F := GF(p);
      R<x> := PolynomialRing(F);
      I := ideal<R | R!fpoly>;
      Q, phi := quo<R|I>;
      elt := phi(Evaluate(R!h,x)*x^p);
      tr := Trace(elt);
      zerolist := [ i : i in [1..#polylist] | Evaluate(R!polylist[i],tr) eq 0 ];
      if #zerolist gt 1 then
        //printf "The prime %o divides the resultant of %o and %o.\n",p,zerolist[1],zerolist[2];
      end if;
      if #zerolist eq 1 then
        //printf "For p = %o, the Frobenius conjugacy class is %o.\n",p,zerolist[1];
        Include(~classesfound,zerolist[1]);
      end if;
      if #zerolist eq 0 then
        printf "For p = %o, we didn't find any classes!\n",p;
      end if;
    end if;
  end if;
end for;

printf "Trying to guess image. Computing maximal subgroups of G.\n";
L := MaximalSubgroups(G);
printf "Done.\n";
printf "Building map from elements to conjugacy classes.\n";
classmap1 := [];
classmap2 := [];
for g in G do
  class := 0;
  for j in [1..#CC] do
    if Index(classlist[j],g) gt 0 then
      class := j;
    end if;
  end for;
  Append(~classmap1,g);
  Append(~classmap2,class);
end for;

for l in [1..#L] do
  //printf "Testing subgroup %o.\n",l;
  classesin := {};
  for g in L[l]`subgroup do
    class := classmap2[Index(classmap1,g)];
    Include(~classesin,class);
  end for;
  if (classesfound subset classesin) then
    printf "The image (might) be contained in maximal subgroup %o.\n",l;
  else
    printf "The image is definitely not contained in maximal subgroup %o.\n",l;
  end if;
end for;