// If an elliptic curve has mod 54 image contained in the subgroup K found
// in the script S3entanglement.txt, then E acquires a 2-torsion point over
// the field that it acquires a cyclic 27-isogeny.

// I'll use Magma's SmallModularCurve databases to finish the computations.
// Magma says that X_0(27) : y^2 + y = x^3 - 7. (Which is, by the way,
// isomorphic to x^3 + y^3 = z^3.) Magma gives the map from X_0(27) to X_0(9)
// that sends (E,C) to (E,3*C) as x9 = phi(x,y) = -3 + y/x + 5/x.

// Magma gives the map from X_0(9) to the j-line as
// j = (x9+9)^3*(x9^3 + 243*x9^2 + 2187*x9 + 6561)^3/x9^9*(x9^2 + 9*x9 + 27)

// First, we want an equation for X_0(27) as a degree 3 cover of Q(X_0(9))

R<y,x,x9> := PolynomialRing(Rationals(),3);
f1 := y^2 + y - (x^3 - 7);
f2 := x*x9 + 3*x - y - 5;
I := ideal<R | f1,f2>;
G := GroebnerBasis(I);

// The second polynomial in the Groebner basis is x^3 - x^2*x9^2 - 6x^2*x9 - 9*x^2 + 9*x*x9 + 27*x - 27 = 0.

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

pol10 := t^3 - t^2*x9^2 - 6*t^2*x9 - 9*t^2 + 9*t*x9 + 27*t - 27;
pol1 := 27*Evaluate(pol10,(1/3)*t + (x9+3)^2/3);
j := (x9+9)^3*(x9^3 + 243*x9^2 + 2187*x9 + 6561)^3/(x9^9*(x9^2 + 9*x9 + 27));
pol2 := t^3 - j*(t+16);

// Now the question is, for which t do these two polynomials define the
// same degree 3 extension?

A1 := Coefficient(pol1,2);
A2 := Coefficient(pol1,1);
A3 := Coefficient(pol1,0);

B1 := Coefficient(pol2,2);
B2 := Coefficient(pol2,1);
B3 := Coefficient(pol2,0);

chk5 := -2*A1*B1;
chk4 := A1^2*B1^2 + 2*A1^2*B2 + 2*A2*B1^2 - 6*A2*B2;
chk3 := -2*A1^3*B1*B2 - 2*A1^3*B3 - 2*A1*A2*B1^3 + 5*A1*A2*B1*B2 + 9*A1*A2*B3 
- 2*A3*B1^3 + 9*A3*B1*B2 - 27*A3*B3;
chk2 := 2*A1^4*B1*B3 + A1^4*B2^2 + 3*A1^2*A2*B1^2*B2 - 9*A1^2*A2*B1*B3 - 
6*A1^2*A2*B2^2 + 2*A1*A3*B1^4 - 9*A1*A3*B1^2*B2 + 27*A1*A3*B1*B3 + A2^2*B1^4 - 
6*A2^2*B1^2*B2 + 9*A2^2*B2^2;
chk1 := -2*A1^5*B2*B3 - 2*A1^3*A2*B1^2*B3 - A1^3*A2*B1*B2^2 + 15*A1^3*A2*B2*B3
- 2*A1^2*A3*B1^3*B2 + 9*A1^2*A3*B1*B2^2 - 27*A1^2*A3*B2*B3 - A1*A2^2*B1^3*B2 + 
9*A1*A2^2*B1^2*B3 + 3*A1*A2^2*B1*B2^2 - 27*A1*A2^2*B2*B3 - 2*A2*A3*B1^5 + 
15*A2*A3*B1^3*B2 - 27*A2*A3*B1^2*B3 - 27*A2*A3*B1*B2^2 + 81*A2*A3*B2*B3;
chk0 := A1^6*B3^2 + A1^4*A2*B1*B2*B3 - 9*A1^4*A2*B3^2 - 2*A1^3*A3*B1^3*B3 + 
A1^3*A3*B1^2*B2^2 + 9*A1^3*A3*B1*B2*B3 - 4*A1^3*A3*B2^3 + A1^2*A2^2*B1^3*B3 - 
9*A1^2*A2^2*B1*B2*B3 + A1^2*A2^2*B2^3 + 27*A1^2*A2^2*B3^2 + A1*A2*A3*B1^4*B2 + 
9*A1*A2*A3*B1^3*B3 - 9*A1*A2*A3*B1^2*B2^2 - 27*A1*A2*A3*B1*B2*B3 + 
18*A1*A2*A3*B2^3 - 4*A2^3*B1^3*B3 + A2^3*B1^2*B2^2 + 18*A2^3*B1*B2*B3 - 
4*A2^3*B2^3 - 27*A2^3*B3^2 + A3^2*B1^6 - 9*A3^2*B1^4*B2 + 27*A3^2*B1^2*B2^2 - 
27*A3^2*B2^3;

S<t> := PolynomialRing(F);

pol := t^6 + chk5*t^5 + chk4*t^4 + chk3*t^3 + chk2*t^2 + chk1*t + chk0;

// If the discriminant of pol is zero, then the degree 3 fields are isomorphic
// if and only if pol has a rational root.

// disc = 0 if x9 = -9. Also, the denominator of disc is 0 if x9 = 0.

// This is a pretty complicated polynomial. Let's use Mark van Hoeij's method
// for finding another degree 6 polynomial that defines the same
// function field.

K := FunctionField(pol);
n := Degree(K);

// Step 1 - Find a basis for OK
OK := MaximalOrderFinite(K);
basis1 := [ OK.i : i in [1..n]];
basis1mat := Matrix(F,[ Eltseq(K!OK.i) : i in [1..n]]);
// Step 2 - Find a basis for Oinfinity
OK2 := MaximalOrderInfinite(K);
basis2 := [ OK2.i : i in [1..Degree(K)]];
basis2mat := Matrix(F,[ Eltseq(K!OK2.i) : i in [1..n]]);
// Step 3 - Write each element of OK as an LC of elements in OK2.
changemat := basis1mat*basis2mat^(-1);
// Step 4 - Let D be the LCM of the denominators of the entries in changemat.
D := LCM([ Denominator(d) : d in Eltseq(changemat)]);
alist := [[Numerator(D*changemat[i][j]) : j in [1..n]] : i in [1..n]];


curb := basis1;
done := false;
while done eq false do
  // Step 5 - Compute max degrees and vectors
  vmat := ZeroMatrix(Rationals(),n,n);
  mlist := [];
  dlist := [];
  for i in [1..Degree(K)] do
    mi := Max([ Degree(a) : a in alist[i]]);
    for j in [1..Degree(K)] do
      vmat[i][j] := Coefficient(alist[i][j],mi);
    end for;
    di := mi - Degree(D);
    Append(~dlist,di);
    Append(~mlist,mi);
  end for;
  // Step 6 - Test linear independence.
  if Rank(vmat) eq n then
    done := true;
  else
    c := Basis(NullSpace(vmat))[1];
    // Step 7 - Update
    maxdi := Max([ dlist[i] : i in [1..#dlist] | c[i] ne 0]);
    ind := 0;
    for i in [1..n] do
      if dlist[i] eq maxdi then
        if c[i] ne 0 then
	  ind := i;
	end if;
      end if;
    end for;
    newb := &+[ c[k]*x9^(dlist[ind]-dlist[k])*curb[k] : k in [1..n]];
    curb[ind] := newb;    
    for j in [1..n] do
      newa := &+[ c[k]*x9^(dlist[ind]-dlist[k])*alist[k][j] : k in [1..n]];
      alist[ind][j] := newa;
    end for;
  end if;
end while;

minpoly := MinimalPolynomial(curb[6]);
minpoly2 := Evaluate(minpoly,t-Coefficient(minpoly,5)/6);
minpoly3 := 458795923^6*Evaluate(minpoly2,108*t/458795923)/108^6;

// So the equation of the genus 4 curve is
// t^6 + (-2*x9^3 - 18x9^2 - 54x9)*t^3 + (x9^6 + 18*x9^5 + 135*x9^4 + 513*x9^3 + 972*x9^2 + 729*x9).

// The map t -> t^3 maps to the following curve.

R<y,x9> := PolynomialRing(Rationals(),2);
pol := y^2 + (-2*x9^3 - 18*x9^2 - 54*x9)*y + (x9^6 + 18*x9^5 + 135*x9^4 + 513*x9^3 + 972*x9^2 + 729*x9);
Y := ProjectiveClosure(Curve(AffineSpace(Rationals(),2),pol));
pt := Y![0,0,1];

E, phi := EllipticCurve(Y,pt);
E2, phi2 := MinimalModel(E);

// So E2 is isomorphic to y^2 = x^3 + 1. Thus, E2 has exactly 6 rational points.

PointSearch(Y,1000);

// The 6 rational points on Y are (-324 : -9 : 1), (0 : 0 : 1), (1 : 0 : 0),
// (-162 : -9 : 1), (0 : -3 : 1) and (-54 : -3 : 1).

FF := FunctionField(Y);
newx9 := FF.2;
jfunc := (newx9+9)^3*(newx9^3 + 243*newx9^2 + 2187*newx9 + 6561)^3/(newx9^9*(newx9^2 + 9*newx9 + 27));

Evaluate(jfunc,Y![-324,-9,1]);
Evaluate(jfunc,Y![0,0,1]);
Evaluate(jfunc,Y![1,0,0]);
Evaluate(jfunc,Y![-162,-9,1]);
Evaluate(jfunc,Y![0,-3,1]);
Evaluate(jfunc,Y![-54,-3,1]);

// The first and fourth points map to j = 0. The second and third points are cusps.
// The fifth and sixth points map to j = -12288000.