
// This is one of two MAGMA files accompanying the paper:

// Title :  Explicit moduli spaces for congruences of elliptic curves
// Author : Tom Fisher
// Date :   April 2018

// In this file "congr-checkpaper.m" we go through all the formulae
// in the paper, in roughly the order of the paper, checking they
// are correct. At the end we also use the formulae to construct
// some examples of pairs of N-congruent elliptic curves over Q.
// Comparing the traces of Frobenius of these curves then provides
// some check on our calculations. We start by loading in the formulae
// from the other file: "congr-formulae.m".

load "congr-formulae.m";
Z := ModularDiagonalQuotientSurface;
R<x,y> := PolynomialRing(K,2);

// Weierstrass equations in Theorem 1.1

W65 := y^2 + 3*T*(T - 2)*x*y + 2*(T - 1)*(T + 2)^2*(T^3 - 2)*y
       - (x^3 - 6*(T - 1)*(T^3 - 2)*x^2);

W73 := y^2 - (x^3 + (4*T^4 + 4*T^3 - 51*T^2 - 2*T - 50)*x^2
       + (6*T + 25)*(52*T^2 - 4*T + 25)*x);

W83 := y^2 - (x^3 - (3*T^2 - 7)*x^2 - 4*T^2*(4*T^4 - 15)*x
       + 4*T^2*(53*T^4 + 81*T^2 + 162));

W85 := y^2 - (x^3 - 2*(T^2 + 19)*x^2
       - (4*T^2 - 49)*(T^4 - 6*T^2 + 25)*x);

W91 := y^2 + (6*T^2 + 3*T + 2)*x*y + T^2*(T + 1)*(4*T^3 + 9*T
       + 9)*y - (x^3 - (16*T^4 + 12*T^3 + 9*T^2 + 6*T + 1)*x^2);

W121 := y^2 + 2*(5*T^2 + 9)*x*y + 96*(T^2 + 3)*(T^2 + 1)^2*y
       - (x^3 + (T^2 + 3)*(11*T^2 + 1)*x^2);

// Weierstrass equations in Theorem 1.2

W87 := y^2 - (x^3 + 2*(4*T^6 - 15*T^4 + 14*T^2 - 1)*x^2
       + (T^2 - 1)^4*(16*T^4 - 24*T^2 + 1)*x);

W92 := y^2 + 3*(4*T^3 + T^2 - 2)*x*y + (T - 1)^3*(T^3 - 1)*(4*T^3
   - 3*T - 7)*y - (x^3 - 3*(T + 1)*(T^3 - 1)*(9*T^2 + 2*T + 1)*x^2);

W101 := y^2 - (3*T - 2)*(6*T^2 - 5*T - 2)*x*y
        - 4*T^2*(T - 1)^2*(4*T^2 - 2*T - 1)*(27*T^3 - 54*T^2 + 16*T + 12)*y
	- (x^3 + T^2*(T - 1)*(27*T^3 - 54*T^2 + 16*T + 12)*x^2);

W103 := y^2 + (T^3 - 8*T^2 - 9*T - 8)*x*y
        + 2*T^2*(T^3 - T^2 - 3*T - 3)*(7*T^2 + 2*T + 3)*y
	- (x^3 + 2*(3*T + 2)*(T^3 - T^2 - 3*T - 3)*x^2);

W111 := y^2 + (T^3 + T)*x*y - (x^3 - (4*T^5 - 17*T^4 + 30*T^3
       - 18*T^2 + 4)*x^2 + T^2*(2*T - 1)*(3*T^2 - 7*T + 5)^2*x);

// Checking the above Weierstrass equation agree with those recorded
// in "congr-formulae.m".

for Neps in [[6,5],[7,3],[8,3],[8,5],[8,7],[9,1],[9,2],
             [10,1],[10,3],[11,1],[12,1]] do
  N,eps := Explode(Neps);
  W := eval Sprintf("W%o%o",N,eps);
  a1,a2,a3,a4,a6 := Explode([MC(W,m): m in [x*y,x^2,y,x,1]]);
  a2 *:= -1;
  a4 *:= -1;
  a6 *:= -1;
  assert W eq (y^2 + a1*x*y + a3*y) - (x^3 + a2*x^2 + a4*x + a6);
  assert EllipticCurve([a1,a2,a3,a4,a6]) eq Z(N,eps);
end for;

// Tables 1 to 4 (as in the paper)

Table1 := [
<[6,5],{*<"I2",2>,<"I3",1>,<"I3",3>,<"I4",1>,<"I4",1>*},1,2,2,20>,
<[7,3],{*<"I1",1>,<"I2",1>,<"I2",2>,<"I2",2>,<"I3",1>,<"I10",1>*},2,2,2,20>,
<[8,3],{*<"I1",2>,<"I2",1>,<"I2",2>,<"I2",2>,<"I3",2>,<"I0*",1>*},1,4,5,20>,
<[8,5],{*<"I2",1>,<"I2",1>,<"I2",2>,<"I2",2>,<"I3",2>,<"I0*",1>*},2,2,4,20>,
<[9,1],{*<"I1",3>,<"I2",1>,<"I2",3>,<"I3",1>,<"I4",1>,<"I0*",1>*},1,3,4,19>,
<[12,1],{*<"I1",8>,<"I4",2>,<"I4",2>*},1,3,5,19>];

Table2 := [
<[8,7],{*<"I2",1>,<"I2",2>,<"I2",2>,<"I3",2>,<"I4",1>,<"I8",1>,<"I8",1>*},          2,1,2,30>,
<[9,2],{*<"I2",3>,<"I3",2>,<"I3",3>,<"I9",1>,<"I0*",1>*},                           1,2,2,29>,
<[10,1],{*<"I2",2>,<"I2",2>,<"I3",3>,<"I5",1>,<"I10",1>,<"IV",1>*},                 1,1,1,28>,
<[10,3],{*<"I1",3>,<"I2",1>,<"I2",2>,<"I2",2>,<"I3",3>,<"I4",1>,<"I4",1>,<"I6",1>*},1,3,4,28>,
<[11,1],{*<"I1",3>,<"I2",1>,<"I2",3>,<"I3",1>,<"I4",1>,<"I4",2>,<"I10",1>*},        2,2,2,28>];

Table3 := [
<[6,5],[0, 2*T^4 - 4*T]>,
<[7,3],[4*T^2 + 20*T + 25, 6*T + 25]>,
<[8,3],[-7, -T^2 + 9, -4*T^2 - 6*T, (4*T^5 - 2*T^4 + 10*T^3 + 6*T^2 + 18*T)/(T - 1)^2]>,
<[8,5],[-4*T^2 + 49, 2*T^3 + 19*T^2 + 60*T + 63]>,
<[9,1],[0, 4*T^4 + 2*T^3 - 2*T^2 , 4*T^4 + 4*T^3 + 9*T^2 + 18*T + 9]>,
<[12,1],[0, -12*T^4 - 24*T^2 - 12, 4*T^6 + 12*T^4 - 4*T^2 - 12]>,
<[8,7],[4*T^6 + 4*T^5 - 9*T^4 - 10*T^3 + 4*T^2 + 6*T + 1]>,
<[9,2],[0, 2*T^5 - 8*T^3 + 4*T^2 + 6*T - 4]>,
<[10,1],[0]>,
<[10,3],[0, 2*T^5 - 4*T^4 - 4*T^3 + 6*T, 4*T^5 - 2*T^4 - 14*T^3 - 18*T^2 - 6*T]>,
<[11,1],[T^4 + 4*T^2 + 4, 3*T^2 - 7*T + 5]>];

Table4 := [
<[8,3],-2,-2*T^4 - 5*T^2 - 9>,
<[8,5],-3,-2*T^3 + T^2 + 18*T - 35>,
<[8,5],-1,16*T^2 - 196>,
<[9,1],-3,-(19/3)*T^4 - 15*T^3 - 9*T^2>,
<[12,1],-3,-12*T^4 - 40*T^2 - 12>,
<[12,1],-1,-16*T^4 - 64*T^2 - 48>,
<[8,7],-3,-4*T^6 - 20*T^5 - 39*T^4 - 36*T^3 - 14*T^2 + 1>,
<[10,3],-3,-7*T^6 - 23*T^5 - 30*T^4 - 15*T^3 - 9*T^2>];

print "Checking Tables 1 and 2";

for dd in Table1 cat Table2 do
  N,eps := Explode(dd[1]);
  E := Z(N,eps);
  LI := LocalInformation(E);
  assert dd[2] eq {* <Sprint(x[5]),Degree(x[1])> :  x in LI *};
  assert dd[3] eq #TorsionSubgroup(E);
  assert exists(xx){x[2] : x in Table3 | x[1] eq [N,eps]};
  MWgens := [Points(E,x)[1]: x in xx];
  reg := Regulator(MWgens);
  printf "(N,eps) = (%o,%o) : Regulator = %o\n",N,eps,reg;
  assert reg ne 0;
  assert dd[4] eq #MWgens; // lower bound for rank over Q(T)
  rk := #MWgens;  
  for xx in Table4 do
    if xx[1] eq [N,eps] then
      L<u> := ext<K|Polynomial(K,[-xx[2],0,1])>;
      P := Points(BaseChange(E,L),xx[3])[1];
      assert Height(P) ne 0;
      rk +:= 1;
    end if;  
  end for;
  assert dd[5] eq rk; // lower bound for rank over Qbar(T)
  rho := rk + 2;
  for x in LI do
    pl,val,f,c,kod := Explode(x);
    m := val + 1 - f; // Ogg's formula
    rho +:= Degree(pl)*(m - 1); // Shioda-Tate formula
  end for;
  assert dd[6] eq rho; // lower bound for geometric Picard number
end for;

print "Section 2 : Proofs via equations for X_E(N,eps)";

print "Case (N,eps) = (7,3)";

R<a,b,x,y,z> := PolynomialRing(Rationals(),5);
F := Evaluate(XE_Equations(7,3,[a,b])[1],[x,y,z]);
S<la> := PolynomialRing(R);
assert Evaluate(F,[la^2*a,la^3*b,la*x,y,la^2*z]) eq la^8*F;
R<x1,x2> := PolynomialRing(K,2);
RR<y> := PolynomialRing(R);
c := (T^2 + 1)*(3*T^2 - 8*T + 3);
h := T^3*(19*T - 42)*x1^2 + 2*T*(T - 3)^2*x1*x2 - 6*(T^2 - 1)*x2^2;
f := 36*T^2*x1^2*x2^2 - 2*(T - 6)*x1*x2^3 + x2^4;
assert Evaluate(F,[y,x1*y,x2,T,y]) eq y^2*(c*y^2 + h*y - f);
quartic := h^2 + 4*c*f;
R<x> := PolynomialRing(K);
C := HyperellipticCurve(Evaluate(quartic,[1,x]));
E := EllipticCurve(C,C![0,T^3*(19*T - 42)]);
TT := (6*T - 3)/(4*T + 4);
E := EllipticCurve([Evaluate(a,TT): a in Eltseq(E)]);
assert IsIsomorphic(Z(7,3),E);

print "Case (N,eps) = (11,1)";

R<a,b,v,w,x,y,z> := PolynomialRing(Rationals(),7);
F := Evaluate(XE_Equations(11,1,[a,b])[1],[v,w,x,y,z]);
S<la> := PolynomialRing(R);
assert Evaluate(F,[la^2*a,la^3*b,la^2*v,la^2*w,la*x,la*y,z]) eq la^6*F;
FF := [ v*z + 2*w*z + x^2 - x*y - y^2,
    a*x*z + b*z^2 - v*x + v*y - 2*w*x,
    a^2*z^2 + 2*a*w*z - 4*a*x^2 - 12*b*x*z - 3*w^2,
    a^2*z^2 + 2*a*v*z - 2*a*x^2 + 4*a*x*y + a*y^2
        - 6*b*x*z + 6*b*y*z + 3*v^2,
    2*a^2*y*z - a*b*z^2 - 4*a*v*x - 2*a*v*y - 6*b*v*z
        - 3*b*x^2 - 12*b*x*y - 3*b*y^2 ];
KK<a,b> := FunctionField(Rationals(),2);
R<v,w,x,y,z> := PolynomialRing(KK,5);
F := Evaluate(F,[a,b,v,w,x,y,z]);
FF1 := [Deriv(F,i): i in [1..5]];
FF2 := [Evaluate(f,[a,b,v,w,x,y,z]): f in FF];
mons := MD(R,2);
mat1 := Matrix(5,#mons,[MC(f,mon): mon in mons,f in FF1]);
mat2 := Matrix(5,#mons,[MC(f,mon): mon in mons,f in FF2]);
M := Solution(mat2,mat1);
assert Determinant(M) eq (4/3)*(4*a^3 + 27*b^2);
assert LCM([Denominator(x): x in Eltseq(M)]) eq 1;
M := Matrix(5,5,[Deriv(f,i): i in [3..7],f in FF]);
R6<U,x1,x2,x3,x4,x5> := PolynomialRing(Rationals(),6);
a := (4*U + 3*x3)*x4 - 3*x5^2;
b := x2*(x1 + x3)*x4 - (4*U + 3*x3)*x4*x5 + 2*x5^3;
xx := [x2*x4 + x4*x5 + x5^2,x3*x4 - x5^2,x5,x4,1];
h := ExactQuotient(4*a^3 + 27*b^2,x4);
M := Matrix(R6,5,5,[Evaluate(M[i,j],[a,b] cat xx): i,j in [1..5]]);
M[2] +:= x5*M[1];
for i in [2..5] do
  M[i] := Vector([ExactQuotient(M[i,j],x4): j in [1..5]]);
end for;  
I := ideal<R6|Minors(M,4)>;
print "Computing the colon ideal (takes a few seconds)";
time J := ColonIdeal(I,ideal<R6|x2*h>);
J1 := EliminationIdeal(J,{U,x1,x2,x3,x4});
R5<U,x1,x2,x3,x4> := PolynomialRing(Rationals(),5);
J1 := ideal<R5|[Evaluate(f,[U,x1,x2,x3,x4,0]): f in Basis(J1)]>;

qq := [ 8*U^2*x2*x4 - 8*U*x1*x2*x3 - 8*U*x1*x3^2 - 8*U*x1*x3*x4 -
    8*U*x2*x3^2 + 12*U*x2*x3*x4 - 8*U*x3^3 - 8*U*x3^2*x4 - 3*x1^2*x2^2
    - 8*x1^2*x2*x3 + 3*x1*x2^3 - 4*x1*x2^2*x3 - 21*x1*x2*x3^2 +
    3*x1*x2*x3*x4 - 4*x1*x3^3 - 4*x1*x3^2*x4 + 3*x2^3*x3 - x2^2*x3^2 -
    13*x2*x3^3 + 9*x2*x3^2*x4 - 4*x3^4 - 3*x3^3*x4, 8*U*x1*x2^2 +
    16*U*x1*x2*x3 - 12*U*x2^2*x4 + 8*U*x2*x3^2 - 8*U*x2*x3*x4 +
    12*x1^2*x2^2 - 6*x1*x2^3 + 26*x1*x2^2*x3 - 6*x1*x2^2*x4 +
    8*x1*x2*x3^2 + 8*x1*x2*x3*x4 + 3*x2^4 - 4*x2^3*x3 + 9*x2^2*x3^2 -
    15*x2^2*x3*x4 + 4*x2*x3^3 + 8*x2*x3^2*x4 + 4*x3^3*x4, 32*U^2*x2*x4
    + 16*U*x1^2*x2 + 32*U*x1^2*x3 + 32*U*x1^2*x4 + 16*U*x1*x2^2 +
    16*U*x1*x2*x3 - 48*U*x1*x2*x4 + 16*U*x1*x3^2 - 16*U*x1*x4^2 -
    16*U*x2^2*x4 - 16*U*x2*x3^2 - 8*U*x2*x3*x4 + 12*U*x2*x4^2 -
    16*U*x3^3 - 32*U*x3^2*x4 - 8*U*x3*x4^2 + 32*x1^3*x2 - 8*x1^2*x2^2
    + 60*x1^2*x2*x3 - 36*x1^2*x2*x4 + 16*x1^2*x3^2 + 16*x1^2*x3*x4 +
    6*x1*x2^3 + 6*x1*x2^2*x4 + 10*x1*x2*x3^2 - 68*x1*x2*x3*x4 +
    6*x1*x2*x4^2 + 8*x1*x3^3 - 8*x1*x3^2*x4 - 8*x1*x3*x4^2 + 6*x2^4 +
    10*x2^3*x3 - 3*x2^3*x4 - 2*x2^2*x3^2 - 8*x2^2*x3*x4 - 26*x2*x3^3 +
    x2*x3^2*x4 + 15*x2*x3*x4^2 - 8*x3^4 - 10*x3^3*x4 - 2*x3^2*x4^2 ];

assert J1 eq ideal<R5|qq>;
XX := Scheme(Proj(R5),J1);
assert IsIrreducible(XX);
assert Dimension(XX) eq 2;
assert Degree(XX) eq 12;

MobiusCoeffs := func<x1,x2,x3,x4|
        4*(2*x1 - x2 + x3), x1*x2 - x2^2 + x3*x4,
	2*(2*x4 - x2), 2*(x1*x2 - x2^2 + x3*x4)>;
QuadForms := func<T,x1,x2,x3,x4| 
       [4*(T - 2)*x1*x2 + 8*x1*x3 + (T - 2)^2*x2^2 + 4*(T - 2)*x2*x3
     + 2*(T^2 - T + 1)*x2*x4 + 4*x3^2 - 4*T*x3*x4,
       8*x1^2 + 16*x1*x3 - 4*(2*T - 1)*x1*x4 - T*(T - 2)*x2^2
     - T^2*x2*x3 - 2*(T^2 - T + 1)*x2*x4 - 2*(T - 4)*x3^2
     - 2*(T^2 + 4*T - 1)*x3*x4]>;
a,b,c,d := MobiusCoeffs(x1,x2,x3,x4);
QQ := QuadForms((a*U + b)/(c*U + d),x1,x2,x3,x4);
assert forall{Q : Q in QQ | Numerator(Q) in J1};

R4<x1,x2,x3,x4> := PolynomialRing(K,4);
QQ := QuadForms(T,x1,x2,x3,x4);
C := Curve(Proj(R4),QQ);
KC<x1,x2,x4> := FunctionField(C);
x3 := 1;
assert forall{Q : Q in QQ | Evaluate(Q,[x1,x2,x3,x4]) eq 0};
a,b,c,d := MobiusCoeffs(x1,x2,x3,x4);
U := (d*T - b)/(-c*T + a);	 
assert forall{q : q in qq | Evaluate(q,[U,x1,x2,x3,x4]) eq 0};
E := EllipticCurve(C,C![0,0,0,1]);
assert IsIsomorphic(Z(11,1),E);

print "Case (N,eps) = (12,1)";

R<a,b,x,y,z,u0,u1,u2> := PolynomialRing(Rationals(),8);
FF := [Evaluate(F,[x,y,z,u0,u1,u2]): F in XE_Equations(12,1,[a,b])];
S<la> := PolynomialRing(R);
subst := [la^2*a,la^3*b,la*x,y,la*z,la^2*u0,la*u1,u2];
assert forall{ i : i in [1..4] |
           Evaluate(FF[i],subst) eq la^[3,2,3,4][i]*FF[i] };    
R<w,x,y> := PolynomialRing(K,3);
v := 2*(w - 2*y)*y + (T + 1)/(T - 1)*(x^2*(y + 1) - (w + 2*y)^2);
xx := [x^2 - 6*y,1,y,v*x,w*x,x];
a := 2*v + w^2 - x^2;
b := (1/6)*(x^2*y - a)*(x^2 - 12*y) + 8*y^3;
FF := [Evaluate(F,[a,b] cat xx): F in FF];
assert [FF[1],FF[2]] eq [0,0];
f := -2^6*(T + 3)^2/(3^4*(T + 1)^2*(T - 1)^6);
g := x^6*(y + 1) - 9*y^2*(x^2 + 4)^2;
h := (T + 1)^2*x^2*y^2 + (T + 2)*(T^2 + 3)^2*x^2
     - 4*(T - 1)*(T + 3)^2*x*y + 4*(T + 3)^2*y^2
     + 12*T*(T + 1)^2*(T + 3)^2;
yy := (144*(T + 1)*y + (T + 3)^2*((T - 3)*x^2
                           + 12*(T + 1)))/(8*(T + 3)*x);
hh := Evaluate(h,[w,x,yy]);
assert Resultant(FF[3],FF[4],w) eq f*x^14*y^2*g^2*hh;
R<x> := PolynomialRing(K);
RR<y> := PolynomialRing(R);
h := Evaluate(h,[0,2*(T + 3)/(T^2 + 3)*x,y]);
F := -(T + 2)*x^4 - (4*T^3 + 5*T^2 + 6*T + 9)*x^2 - 3*T*(T^2 + 3)^2;
assert IsSquare(Discriminant(h)/F);

R<x1,x2,x3> := PolynomialRing(K,3);
G := x1^2*x2 + (T^2 + 2)*x1^2*x3 + 2*T*x1*x2^2 - 2*T*x1*x2*x3
            + T^2*x2^3 + 3*x2^2*x3 + 3*T^2*x2*x3^2 + 9*x3^3;
G0 := x1^2*x2 + 2*x1^2*x3 + 3*x2^2*x3 + 9*x3^3;
G1 := 2*x1*x2*(x2 - x3);
G2 := x1^2*x3 + x2^3 + 3*x2*x3^2;
assert G eq G2*T^2 + G1*T + G0;
assert Evaluate(G1^2 - 4*G0*G2,[x,T,1]) eq 4*F;
C := Curve(Proj(R),G);
E := EllipticCurve(C,C![1,0,0]);
assert IsIsomorphic(Z(12,1),E);

print "Section 3 : Degree 3 covers of K3 surfaces";

// The equations for Z^*(N,eps) in Theorem 3.1

function CuspidalCubics(N,eps,uu)
  u,v,w := Explode(uu);
  case [N,eps] :
    when [3,2] : 
      Fplus := u*(u + 3*v + w)^2 + 4*v^3;
      Fminus := u*(u + 3*v - w)^2 + 4*v^3;
    when [5,1] :
      Fplus := u*(u^2 - 11*u*v - v^2) + w^2*(12*u + v)
                     + 2*w*(3*u^2 - 4*u*v + 4*w^2 );
      Fminus := u*(u^2 - 11*u*v - v^2) + w^2*(12*u + v)
                     - 2*w*(3*u^2 - 4*u*v + 4*w^2 );
    when [5,2] :   
      Fplus := u^2*(11*v + 8*w) + w^2*(8*u - v + 4*w)
                     + 2*u*(2*v - w)*(4*u - v + 4*w);
      Fminus := u^2*(11*v + 8*w) + w^2*(8*u - v + 4*w)
                     - 2*u*(2*v - w)*(4*u - v + 4*w);
  end case;
  return Fplus,Fminus;
end function;

R<x,a,b> := PolynomialRing(Rationals(),[1,2,3]);
Delta := -4*a^3 - 27*b^2;
f := x^3 + a*x + b;
fx := 3*x^2 + a;
g := 3*a*x^4 + 18*b*x^3 - 6*a^2*x^2 - 6*a*b*x - a^3 - 9*b^2;
h := 3*a*x^2 + 9*b*x - a^2;
j := 27*b*x^3 - 18*a^2*x^2 - 27*a*b*x - 2*a^3 - 27*b^2;
k := x^9 + 12*a*x^7 + 84*b*x^6 - 42*a^2*x^5 - 84*a*b*x^4
   - 20*(a^3 + 12*b^2)*x^3 + 12*a^2*b*x^2
   - 3*a*(5*a^3 + 32*b^2)*x - 4*b*(3*a^3 + 20*b^2);

assert j^2 eq -4*h^3 - 27*Delta*f^2;
assert k eq f^3 + f*j + 4*Delta*f + 3*g*(x*fx - 2*f);

print "Case (N,eps) = (3,2)";

S<xi,eta> := PolynomialRing(R,2);
D := -27*a*xi^4 - 54*b*xi^3*eta - 18*a^2*xi^2*eta^2
                         - 54*a*b*xi*eta^3 + (a^3 - 27*b^2)*eta^4;
D1 := HessePolynomials(3,2,[-48*a,-864*b]:Variables := [12*xi,eta]);
assert D eq (1/(2^12*3^3))*D1;
Hmat := Matrix(2,2,[Deriv(Deriv(D,i),j): i,j in [1..2]]);
A := (1/108)*Determinant(Hmat);
Jmat := Matrix(2,2,[Deriv([D,A][i],j): i,j in [1..2]]);
B := (1/36)*Determinant(Jmat);
assert -4*A^3 - 27*B^2 eq 16*(4*a^3 + 27*b^2)^2*D^3;
D := Evaluate(D,[x,1]);
A := Evaluate(A,[x,1]);
B := Evaluate(B,[x,1]);
assert D eq fx^3 - 27*f^2;
assert D eq j - 3*fx*h;
subst := [D,fx*h,h^3,3^6*Delta*f^4];
RR<u,v,r,s> := PolynomialRing(Rationals(),[1,1,2,3]);
F := (4*r + (u + 3*v)^2)*(r*u - v^3) - r*s;
assert IsHomogeneous(F);
assert Evaluate(F,subst) eq 0;

// Formulae for the inverse map (Remark 3.2)
invmap32 := func<u,v,r|
  r + v^2,
  -3*r*(r + u*v + 2*v^2),
  r*(u + 3*v)*(r*u + v^3) + 2*r^2*(r + 3*v^2) >;
xx,aa,bb := invmap32(u,v,r);
uu := [Evaluate(s,[xx,aa,bb]): s in subst];
la := -27*r*(r*u - v^3)^3;
vv := [la^[1,1,2,3][i]*[u,v,r,s][i] : i in [1..4]];
assert forall{i : i in [1..4] | uu[i] - vv[i] in ideal<RR|F>};
u,v,r,s := Explode(subst);
xx,aa,bb := invmap32(u,v,r);
la := 9*f*h^2;
assert [xx,aa,bb] eq [la*x,la^2*a,la^3*b];

// The involution that switches the two elliptic curves
// (We replace r by the other root of the quadratic)
r1 := -(u*(u + 3*v)^2 - 4*v^3 - s)/(4*u) - r;
_,aa,bb := invmap32(u,v,r1);
mu := 9/4*fx^2*f*j/D;
assert [aa,bb] eq [mu^2*A,mu^3*B];
assert la*mu eq (9/2*fx*f*h)^2*(u + 3*v)/u;

// Proof of Theorem 3.1
RR<u,v,w> := PolynomialRing(Rationals(),3);
Fplus,Fminus := CuspidalCubics(3,2,[u,v,w]);
P<r> := PolynomialRing(RR);
poly := Evaluate(F,[u,v,r,u*w^2]);
assert Discriminant(poly) eq Fplus*Fminus;

print "Case (N,eps) = (5,1)";

S<la,mu> := PolynomialRing(R,2);
D := la^12 + 22*a*la^10*mu^2 + 220*b*la^9*mu^3 - 165*a^2*la^8*mu^4
  - 528*a*b*la^7*mu^5 - 220*(a^3 + 12*b^2)*la^6*mu^6 + 264*a^2*b*la^5*mu^7
  - 165*a*(5*a^3 + 32*b^2)*la^4*mu^8 - 880*b*(3*a^3 + 20*b^2)*la^3*mu^9
  + 22*a^2*(25*a^3 + 168*b^2)*la^2*mu^10 + 20*(19*a^4*b + 128*a*b^3)*la*mu^11
  + (125*a^6 + 1792*a^3*b^2 + 6400*b^4)*mu^12;
D1 := HessePolynomials(5,1,[-48*a,-864*b]:Variables := [12*la,mu]);
assert D eq (1/(2^24*3^12))*D1;
Hmat := Matrix(2,2,[Deriv(Deriv(D,i),j): i,j in [1..2]]);
A := (1/5808)*Determinant(Hmat);
Jmat := Matrix(2,2,[Deriv([D,A][i],j): i,j in [1..2]]);
B := (1/360)*Determinant(Jmat);
assert 4*A^3 + 27*B^2 eq (4*a^3 + 27*b^2)*D^5;
D := Evaluate(D,[x,1]);
A := Evaluate(A,[x,1]);
B := Evaluate(B,[x,1]);
assert D eq 4*k*f - 3*(f^2 + g)^2 + 32*(f^2 + g)*Delta;
subst := [4*f,2*(f^2 + g),16*Delta,4*k,D];
RR<t,u,v,r,s> := PolynomialRing(Rationals(),[1,2,2,3,4]);
RHS1 := u*(u^2 - 11*u*v - v^2) + (12*u + v)*s;
RHS2 := 3*u^2 - 4*u*v + 4*s;
eqns := [r^2 + s*t^2 - RHS1,r*t - RHS2];
assert forall{ F : F in eqns | IsHomogeneous(F) };
assert forall{ F : F in eqns | Evaluate(F,subst) eq 0 };
RHS := RHS1^2 - 4*s*RHS2^2;

// Formulae for the inverse map (Remark 3.2)
invmap51 := func<u,v,rt,tt|
  32*u - v + 5*tt,
  -3*(8*u - v)*(32*u - v) - 288*rt + 30*(28*u + v)*tt - 75*tt^2,
  -2*(32*u - v)^2*(4*u + v) - 144*(32*u - v)*rt
       + 6*(32*u - v)*(88*u - 5*v)*tt + 1008*rt*tt
       - 150*(28*u + v)*tt^2 + 250*tt^3 >;
xx,aa,bb := invmap51(u,v,r*t,t^2);
uu := [Evaluate(s,[xx,aa,bb]): s in subst];
la := 1728*(-(32*u - v)*r + (125*u^2 - 4*s)*t);
vv := [la^[1,2,2,3,4][i]*[t,u,v,r,s][i] : i in [1..5]];
assert forall{i : i in [1..5] | uu[i] - vv[i] in ideal<RR|eqns>};
t,u,v,r,s := Explode(subst);
xx,aa,bb := invmap51(u,v,r*t,t^2);
la := 24*Derivative(f^2 + g,x);
assert [xx,aa,bb] eq [la*x,la^2*a,la^3*b];

// The involution that switches the two elliptic curves
// (We replace r and t by sqrt(s)*t and 1/sqrt(s)*r)
_,aa,bb := invmap51(u,v,r*t,r^2/s);
mu := -2*Derivative(Derivative(k,x),x)/D; 
assert [aa,bb] eq [mu^2*A,mu^3*B];
assert la*mu eq -144*(35*u^2 - 5*u*v + 4*s)/s;

// Proof of Theorem 3.1
RR<u,v,w> := PolynomialRing(Rationals(),3);
Fplus,Fminus := CuspidalCubics(5,1,[u,v,w]);
assert Evaluate(RHS,[0,u,v,0,w^2]) eq Fplus*Fminus;

print "Case (N,eps) = (5,2)";

S<la,mu> := PolynomialRing(R,2);
D := (125*a^3 - 432*b^2)*la^12 + 2430*a^2*b*la^11*mu
   - 22*a*(25*a^3 - 378*b^2)*la^10*mu^2
   - 110*b*(11*a^3 - 108*b^2)*la^9*mu^3
   - 165*a^2*(5*a^3 - 27*b^2)*la^8*mu^4
   - 132*a*b*(53*a^3 - 189*b^2)*la^7*mu^5
   + 220*(a^6 - 123*a^3 *b^2 + 81*b^4)*la^6*mu^6
   + 132*a^2*b*(19*a^3 - 297*b^2)*la^5*mu^7
   - 165*(a^7 - 26*a^4*b^2 + 189*a*b^4)*la^4*mu^8
   - 110*(3*a^6*b - 34*a^3*b^3 + 135*b^5)*la^3*mu^9
   - 22*a^2*(a^3 - 3*b^2)*(a^3 + 27*b^2)*la^2*mu^10
   - 10*a*b*(5*a^6 + 82*a^3*b^2 + 189*b^4)*la*mu^11
   + (a^9 - a^6*b^2 - 181*a^3*b^4 - 675*b^6)*mu^12;
D1 := HessePolynomials(5,2,[-48*a,-864*b]:Variables := [12*la,mu]);
assert D eq (1/(2^36*3^15))*D1;
Hmat := Matrix(2,2,[Deriv(Deriv(D,i),j): i,j in [1..2]]);
A := (1/1452)*Determinant(Hmat);
Jmat := Matrix(2,2,[Deriv([D,A][i],j): i,j in [1..2]]);
B := (-1/180)*Determinant(Jmat);
assert -4*A^3 - 27*B^2 eq 16*(4*a^3 + 27*b^2)^2*D^5;
D := Evaluate(D,[x,1]);
A := Evaluate(A,[x,1]);
B := Evaluate(B,[x,1]);
assert D eq 16*Delta*f^4 - g^3 + 4*(2*g^3 - g^2*j - 4*Delta*f^2*g);
subst := [Delta*f^4,Delta*f^2*g,g^3,2*g^3 - g^2*j - 4*Delta*f^2*g];
RR<r,s,v,w> := PolynomialRing(Rationals(),4);
F := r*(4*s - 2*v + w)^2 + 27*r*s*v + s*w^2 - s^2*(v - 4*w);
assert IsHomogeneous(F);
assert Evaluate(F,subst) eq 0;

// Formulae for the inverse map (Remark 3.2)
invmap52 := func<r,s,v,w|
  4*r*s + 4*r*v + r*w - s^2,
  3*(8*r*s^3 + 4*r*s^2*v + 6*r*s^2*w + r*s*w^2 - s^4),
  r^2*s*(16*s^3 - 8*s^2*v - 24*s^2*w - 40*s*v*w - 15*s*w^2 + 4*v*w^2
     - 2*w^3) + r*s^3*(24*s^2 + 8*s*v + 34*s*w + 7*w^2) - 2*s^6 >;	
xx,aa,bb := invmap52(r,s,v,w);  
uu := [Evaluate(s,[xx,aa,bb]): s in subst];
M := Matrix(RR,2,4,[uu,[r,s,v,w]]);
assert Minors(M,2) subset ideal<RR|F>;
r,s,v,w := Explode(subst);
xx,aa,bb := invmap52(r,s,v,w);  
la := (3/2)*Delta*f^4*g^2*Derivative(g,x);
assert [xx,aa,bb] eq [la*x,la^2*a,la^3*b];
assert D eq 16*r - v + 4*w;

// The involution that switches the two elliptic curves
// (We replace s by the other root of the quadratic)
s1 := -(27*r*v + w^2 - 8*r*(2*v - w))/(16*r - v + 4*w) - s;
_,aa,bb := invmap52(r,s1,v,w);  
mu0 := (150*a^4 - 81*a*b^2)*x^11 + (2745*a^3*b + 486*b^3)*x^10 +
    (-650*a^5 + 15660*a^2*b^2)*x^9 + (-6705*a^4*b + 44955*a*b^3)*x^8 +
    (300*a^6 - 31050*a^3*b^2 + 47385*b^4)*x^7 + (4242*a^5*b -
    55566*a^2*b^3)*x^6 + (-420*a^7 + 11088*a^4*b^2 - 32319*a*b^4)*x^5
    + (-2010*a^6*b + 3780*a^3*b^3 - 18225*b^5)*x^4 + (30*a^8 -
    2445*a^5*b^2 - 2835*a^2*b^4)*x^3 + (-75*a^7*b - 3360*a^4*b^3 -
    11745*a*b^5)*x^2 + (14*a^9 + 156*a^6*b^2 - 819*a^3*b^4 -
    6075*b^6)*x + 11*a^8*b + 169*a^5*b^3 + 594*a^2*b^5;
mu := 3*f^4*g^2*(g*j + 4*Delta*f^2)*mu0/D^2;
assert [aa,bb] eq [mu^2*A,mu^3*B];
assert la*mu eq (3*f*g/D)^2*r*(2*v - w)*(-4*r*v + 3*r*w + w^2);

// Proof of Theorem 3.1
RR<u,v,w> := PolynomialRing(Rationals(),3);
Fplus,Fminus := CuspidalCubics(5,2,[u,v,w]);
P<s> := PolynomialRing(RR);
PP<r> := PolynomialRing(P);
F1 := Evaluate(F,[r,s,v,w]);
F2 := r*(16*r - v + 4*w) - (4*r - u)^2;
assert Discriminant(Resultant(F1,F2)) eq Fplus*Fminus;

// The involution that switches Fplus and Fminus
uu := u*(v - 4*w)/(8*u - (v - 4*w));
assert F2 eq (8*u - v + 4*w)*r - u^2;
assert u^2/(8*u - v + 4*w) eq uu^2/(8*uu - v + 4*w);
assert Evaluate(Fplus,[uu,v,w]) eq (uu/u)^2*Fminus;
assert Evaluate(Fminus,[uu,v,w]) eq (uu/u)^2*Fplus;

print "Checking Remark 3.4";

j3split := 27*(T - 3)^3*(T + 1)^3/T^3;
j5split := (T + 5)^3*(T^2 - 5)^3*(T^2 + 5*T + 10)^3/(T^2 + 5*T + 5)^5;
j5nonsplit := 125*T*(2*T + 1)^3*(2*T^2 + 7*T + 8)^3/(T^2 + T - 1)^5;

// Comparing to formulae for X_0(N^2)

j9 := (T + 3)^3*(T^3 + 9*T^2 + 27*T + 3)^3/(T*(T^2 + 9*T + 27));
j25 := (T^10 + 10*T^9 + 55*T^8 + 200*T^7 + 525*T^6 + 1010*T^5
         + 1425*T^4 + 1400*T^3 + 875*T^2 + 250*T + 5)^3
                /(T*(T^4 + 5*T^3 + 15*T^2 + 25*T + 25));

E := EllipticCurveWithjInvariant(j9);
E1 := IsogenyFromKernel(E,f)
   where f is Factorization(DivisionPolynomial(E,3))[1,1];
assert jInvariant(E1) eq Evaluate(j3split,(1/3)*(T + 27/T) + 3);

E := EllipticCurveWithjInvariant(j25);
E1 := IsogenyFromKernel(E,f)
   where f is Factorization(DivisionPolynomial(E,5))[1,1];
assert jInvariant(E1) eq Evaluate(j5split,T + 5/T);

// Parametrisations of the cuspidal cubics {Fplus = 0}

// Case (N,eps) = (3,2)
param := [ -1, T^2, 2*T^3 - 3*T^2 + 1];
assert CuspidalCubics(3,2,param) eq 0;
param := [Evaluate(f,T/6 - 1/2): f in param];
u,v,w := Explode(param);
PP<r> := PolynomialRing(K);
poly := (4*r + (u + 3*v)^2)*(r*u - v^3) - r*u*w^2;
assert #Roots(poly) eq 1;
r := Roots(poly)[1,1];
_,a,b := invmap32(u,v,r);
assert jInvariant(EllipticCurve([a,b])) eq j3split;

// Case (N,eps) = (5,1)
param := [ T^2*(2*T + 1), 1, -T*(T^2 - T - 1) ];
assert CuspidalCubics(5,1,param) eq 0;
param := [Evaluate(f,1/(T + 2)): f in param];
u,v,w := Explode(param);
cc := -(3*u^2 - 4*u*v + 4*w^2)/w;
u,v,w := Explode([f/cc : f in param]);
r,s,t := Explode([-w,w^2,1]);
assert r^2 + s*t^2 eq u*(u^2 - 11*u*v - v^2) + (12*u + v)*s;
assert r*t eq 3*u^2 - 4*u*v + 4*s;
_,a,b := invmap51(u,v,r*t,t^2);
assert jInvariant(EllipticCurve([a,b])) eq j5split;

// Case (N,eps) = (5,2)
param := [ T + 2, -T^3, 2*T^2 + 3*T ];
assert CuspidalCubics(5,2,param) eq 0;
param := [Evaluate(f,-2*(T + 2)): f in param];
u,v,w := Explode(param);
PP<s> := PolynomialRing(K);
poly := (4*u - v + 4*w)^2*s^2 + (u^2*(11*v + 8*w)
         + (8*u - v + 4*w)*w^2)*s + u^2*(2*v - w)^2;
assert #Roots(poly) eq 1;
s := Roots(poly)[1,1];
r := u^2/(8*u - v + 4*w);
assert r*(4*s - 2*v + w)^2 + 27*r*s*v + s*w^2 - s^2*(v - 4*w) eq 0;
assert r*(16*r - v + 4*w) eq (4*r - u)^2;
_,a,b := invmap52(r,s,v,w);
assert jInvariant(EllipticCurve([a,b])) eq j5nonsplit;

print "Checking the equations for Z(6,5), Z(10,1), Z(10,3)";

// Case (N,eps) = (6,5)
PP<x> := PolynomialRing(K);
subst := [x,1,(T^3 - 1)*x + 3*(T - 1)];
Fplus,Fminus := CuspidalCubics(3,2,subst);
assert exists(f){f[1] : f in Factorization(Fplus) | f[2] eq 2};
C := HyperellipticCurve(ExactQuotient(Fplus,f^2)*Fminus);
assert IsIsomorphic(EllipticCurve(C),Z(6,5));

// Case (N,eps) = (10,1)
PP<x> := PolynomialRing(K);
subst := [x,1,((T - 2)*x - T*(T - 1)^2)/(-2*(T - 1))];
Fplus,Fminus := CuspidalCubics(5,1,subst);
assert exists(f){f[1] : f in Factorization(Fplus) | f[2] eq 2};
C := HyperellipticCurve(ExactQuotient(Fplus,f^2)*Fminus);
assert IsIsomorphic(EllipticCurve(C),Z(10,1));

// Case (N,eps) = (10,3)
PP<x> := PolynomialRing(K);
subst := [x,1,(T^3*x - (T + 1))/T^2];
Fplus,Fminus := CuspidalCubics(5,2,subst);
assert exists(f){f[1] : f in Factorization(Fplus) | f[2] eq 2};
C := HyperellipticCurve(ExactQuotient(Fplus,f^2)*Fminus);
assert IsIsomorphic(EllipticCurve(C),Z(10,3));

print "Proof of Lemma 3.5";

KK<t,a,b> := FunctionField(Rationals(),3);
E1 := EllipticCurve([a,b]);
D,c4,c6 := HessePolynomials(2,1,cInvariants(E1):Variables := [12*t,1]);
E2 := EllipticCurve([-27*c4,-54*c6]);
j1 := jInvariant(E1);
j2 := jInvariant(E2);
m := (1/2^5)*b*c6/(t^3 + a*t + b)/(4*a^3 + 27*b^2);
x := -a*(2*a*t + 3*b)*c4/(t^3 + a*t + b)/(4*a^3 + 27*b^2);
assert (j1 - 1728)*(j2 - 1728) eq m^2;
assert x^3 - 3*j1*j2*x - 2*j1*j2*(m + 1728) eq 0;

KK<j1,j2> := FunctionField(Rationals(),2);
cc := [j1/(j1 - 1728),j1/(j1 - 1728)];
E1 := EllipticCurve([-27*cc[1],-54*cc[2]]);
assert jInvariant(E1) eq j1;
PP<m,x> := PolynomialRing(KK,2);
LL<m,x> := quo<PP|m^2 - (j1 - 1728)*(j2 - 1728),
                  x^3 - 3*j1*j2*x - 2*j1*j2*(m + 1728)>;
t := (x + j1)*(1728 - m - j1)/(12*(j1 - j2)*(j1 - 1728));
D,c4,c6 := HessePolynomials(2,1,cc:Variables := [12*t,1]);
E2 := EllipticCurve([-27*c4,-54*c6]);
assert jInvariant(E2) eq j2;

print "Checking the degree 3 covers are the ones we want";

print "Case (N,eps) = (6,5)";

KK<u,v,w> := FunctionField(Rationals(),3);
Fplus,Fminus := CuspidalCubics(3,2,[u,v,w]);
PP<y> := PolynomialRing(KK);
LL<y> := quo<PP|y^2 - Fplus*Fminus>;
P<r> := PolynomialRing(KK);
F := (4*r + (u + 3*v)^2)*(r*u - v^3) - r*u*w^2;
r1,r2 := Explode([x[1] : x in Roots(F,LL)]);
_,a1,b1 := invmap32(u,v,r1);
_,a2,b2 := invmap32(u,v,r2);
j1 := jInvariant(EllipticCurve([a1,b1]));
j2 := jInvariant(EllipticCurve([a2,b2]));
G4 := u^4 + 4*u^3*v + 6*u^2*v^2 - u^2*w^2 + 4*u*v^3 - 2*u*v*w^2 + v^4;
G6 := 2*u^6 + 12*u^5*v + 30*u^4*v^2 + u^4*w^2 + 40*u^3*v^3 +
    12*u^3*v*w^2 + 2*u^3*w^3 + 30*u^2*v^4 + 21*u^2*v^2*w^2 - u^2*w^4 +
    12*u*v^5 + 10*u*v^3*w^2 + 2*v^6;
m := 2^5*3^3*G6/(u^3*w^3) - 1728;
H := 2^5*3^3*G4/(u^3*w^3);
assert (j1 - 1728)*(j2 - 1728) eq m^2;
assert j1*j2 eq 4*G4*H^2;
assert j1*j2*(m + 1728) eq 4*G6*H^3;
PP<t> := PolynomialRing(KK);
LL<t> := quo<PP|(t^3 - 1)*u + 3*(t - 1)*v - w>;
x0 := 2*u*(u + v)*t^2 + 2*u*(u + v - w)*t + 4*(u + v)*v;
assert x0^3 - 12*G4*x0 - 8*G6 eq 0;

print "Case (N,eps) = (10,1)";

KK<u,v,w> := FunctionField(Rationals(),3);
Fplus,Fminus := CuspidalCubics(5,1,[u,v,w]);
PP<y> := PolynomialRing(KK);
LL<y> := quo<PP|y^2 - Fplus*Fminus>;
RHS1 := u*(u^2 - 11*u*v - v^2) + (12*u + v)*w^2;
rt := 3*u^2 - 4*u*v + 4*w^2;
_,a1,b1 := invmap51(u,v,rt,(RHS1 - y)/(2*w^2));
_,a2,b2 := invmap51(u,v,rt,(RHS1 + y)/(2*w^2));
j1 := jInvariant(EllipticCurve([a1,b1]));
j2 := jInvariant(EllipticCurve([a2,b2]));
G4 := u^4 + 228*u^3*v + 494*u^2*v^2 - 8*u^2*w^2 - 228*u*v^3 -
        312*u*v*w^2 + v^4 + 208*v^2*w^2 + 16*w^4;
G6 := u^6 - 522*u^5*v - 10005*u^4*v^2 - 12*u^4*w^2 + 2484*u^3*v*w^2 -
        10005*u^2*v^4 + 4692*u^2*v^2*w^2 + 48*u^2*w^4 + 522*u*v^5 +
        14076*u*v^3*w^2 - 2448*u*v*w^4 + v^6 - 552*v^4*w^2 -
        4488*v^2*w^4 + 1728*v*w^5 - 64*w^6;
m := G6/(v*w^5) - 1728;
H := G4/(v*w^5);
assert (j1 - 1728)*(j2 - 1728) eq m^2;
assert j1*j2 eq G4*H^2; 
assert j1*j2*(m + 1728) eq G6*H^3; 
PP<t> := PolynomialRing(KK);
LL<t> := quo<PP|(t - 2)*u - t*(t - 1)^2*v + 2*(t - 1)*w>;
x0 := -3*v*(u + 7*v - 2*w)*t^2 + 6*v*(4*u + 3*v - 3*w)*t
 	      + 2*(u^2 + v^2 + 18*v*w - 4*w^2);
assert x0^3 - 3*G4*x0 - 2*G6 eq 0;

print "Case (N,eps) = (10,3)";

KK<u,v,w> := FunctionField(Rationals(),3);
Fplus,Fminus := CuspidalCubics(5,2,[u,v,w]);
PP<y> := PolynomialRing(KK);
LL<y> := quo<PP|y^2 - Fplus*Fminus>;
r := u^2/(8*u - v + 4*w);
PP<s> := PolynomialRing(LL);
F := (4*u - v + 4*w)^2*s^2 + (u^2*(11*v + 8*w)
         + (8*u - v + 4*w)*w^2)*s + u^2*(2*v - w)^2;
s1,s2 := Explode([x[1] : x in Roots(F,LL)]);
_,a1,b1 := invmap52(r,s1,v,w);  
_,a2,b2 := invmap52(r,s2,v,w);  
j1 := jInvariant(EllipticCurve([a1,b1]));
j2 := jInvariant(EllipticCurve([a2,b2]));
G6 := 640*u^4*v^2 - 768*u^4*v*w - 72*u^3*v^3 - 240*u^3*v^2*w +
    320*u^3*v*w^2 + 9*u^2*v^4 - 6*u^2*v^3*w - 160*u^2*v^2*w^2 +
    160*u^2*v*w^3 + 64*u^2*w^4 - 16*u*v*w^4 + 64*u*w^5 + v^2*w^4 -
    8*v*w^5 + 16*w^6;
G9 := 6912*u^7*v^2 - 1376*u^6*v^3 - 14976*u^6*v^2*w + 9828*u^5*v^4 -
    32544*u^5*v^3*w + 49728*u^5*v^2*w^2 - 18432*u^5*v*w^3 -
    2673/2*u^4*v^5 + 10278*u^4*v^4*w - 30264*u^4*v^3*w^2 +
    36960*u^4*v^2*w^3 - 13056*u^4*v*w^4 + 27*u^3*v^6 - 432*u^3*v^5*w +
    3240*u^3*v^4*w^2 - 10224*u^3*v^3*w^3 + 10752*u^3*v^2*w^4 -
    3840*u^3*v*w^5 - 512*u^3*w^6 - 81/2*u^2*v^5*w^2 + 369*u^2*v^4*w^3
    - 1068*u^2*v^3*w^4 + 1200*u^2*v^2*w^5 - 768*u^2*v*w^6 -
    768*u^2*w^7 - 24*u*v^2*w^6 + 192*u*v*w^7 - 384*u*w^8 + v^3*w^6 -
    12*v^2*w^7 + 48*v*w^8 - 64*w^9;
m := 2^6*(8*u - v + 4*w)*G9/(u^3*v^2*(4*u - v + 4*w)^5) - 1728;
H := 2^6*(8*u - v + 4*w)*G6/(u^3*v^2*(4*u - v + 4*w)^5);
assert (j1 - 1728)*(j2 - 1728) eq m^2;
assert j1*j2 eq G6*H^2; 
assert j1*j2*(m + 1728) eq G9*H^3;
PP<t> := PolynomialRing(KK);
LL<t> := quo<PP|t^3*u - (t + 1)*v - t^2*w>;
x0 := 3*u^2*(8*u - 3*v - 4*w)*t^2 + 12*u*(2*u*v - 4*u*w + v*w)*t
          - 16*u^2*v + 6*u*v^2 + 8*u*w^2 - v*w^2 + 4*w^3;
assert x0^3 - 3*G6*x0 - 2*G9 eq 0;

print "Section 4 : Computing the Picard numbers";

// Table 5 (as in the paper, but with a final column recording Delta_p)

P<x> := PolynomialRing(Rationals());
Table5 := [
<[9,1],[
<5,(x - 1)^16*(x + 1)^2*(x^2 + x + 1)*(x^2 + 7/5*x + 1),51>,
<7,(x - 1)^18*(x + 1)^2*(x^2 + 10/7*x + 1),6>]>,
<[12,1],[
<5,(x - 1)^16*(x + 1)^4*(x^2 + 6/5*x + 1),1>,
<11,(x - 1)^12*(x + 1)^8*(x^2 + 6/11*x + 1),7>]>,
<[9,2],[
<7,(x - 1)^24*(x + 1)^2*(x^2 + x + 1)^2*(x^2 + 10/7*x + 1)*(x^2 + 13/7*x + 1),2>,
<13,(x - 1)^24*(x^2 + x + 1)^3*(x^2 + 1/13*x + 1)*(x^2 + 25/13*x + 1),17>]>,
<[10,1],[
<7,(x - 1)^24*(x + 1)^2*(x^2 + x + 1)^2*(x^2 + 10/7*x + 1)^2,1>,
<17,-(x - 1)^25*(x + 1)^5*(x^2 - 2/17*x + 1)*(x^2 + 25/17*x + 1),118>]>,
<[10,3],[
<31,(x - 1)^24*(x + 1)^2*(x^2 + x + 1)^2*(x^2 + 46/31*x + 1)*(x^2 + 58/31*x + 1),10>,
<37,(x - 1)^28*(x + 1)^2*(x^2 + 70/37*x + 1)^2,1>]>,
<[11,1],[
<23,(x - 1)^28*(x + 1)^2*(x^2 + 42/23*x + 1)*(x^2 + 45/23*x + 1),2002>,
<53,(x - 1)^28*(x^2 + x + 1)*(x^2 + 25/53*x + 1)*(x^2 + 70/53*x + 1),1441>]>];

// Extra generators for the Mordell-Weil groups mod p

ExtraGenerators := [
<[9,2],7,[[ 5*T^5 + 6*T^4 + 4*T^2 + 6*T, T^7 + 3*T^6 + 6*T^3 + 2*T^2 +
2*T]]>,
<[9,2],13,[[ 4*T^6 + 8*T^5 + 3*T^4 + 7*T^3 + 5*T^2 + 10*T + 2, 10*T^9
         + 4*T^8 + 5*T^6 + 5*T^5 + 12*T^3 + 4*T^2 + 12 ]]>,
<[10,1],7, [[ 6*T^6 + 6*T^4 + 4*T^3 + 5*T^2, 4*T^9 + 6*T^8 + 6*T^7 +
T^6 + T^5 + 3*T^4], [ T^6 + 5*T^5 + 6*T^4 + 4*T^3 + 5*T^2, 2*T^9 +
6*T^8 + 2*T^7 + T^6 + 3*T^4 ]]>,		    
<[10,1],17,[[ 16*T^6 + 13*T^5 + 6*T^4 + 4*T^3 + 12*T^2, 4*T^9 + 2*T^8
+ 5*T^7 + 8*T^5 + 15*T^4 ], [ (6*T^8 + 8*T^7 + 2*T^6 + 5*T^5 + 8*T^4 +
4*T^3 + T^2)/(T + 6)^2 ]]>,
<[10,3],31,[[ 20*T^4 + 13*T^3 + 30*T^2 + 6*T, 5*T^5 + 4*T^4 + 29*T^3 +
4*T^2 + 5*T ], [ (7*T^6 + 12*T^5 + 9*T^4 + 19*T^2 + 13*T + 4)/(T +
29)^2 ]]>,
<[10,3],37,[[ 36*T^4 + 11*T^3 + 4*T^2 , 26*T^5 + 34*T^4 + 2*T^3 +
15*T^2 ], [ 6*T^4 + 5*T^3 + T^2 + 26*T + 32, 29*T^5 + 35*T^4 + 2*T^3 +
15*T^2 + 19*T + 10 ]]>,
<[11,1],23, [[ 16*T^2 + 5*T + 5, 21*T^3 + 15*T^2 + 3*T + 18 ], [
(18*T^6 + 5*T^5 + 5*T^4 + 22*T^3 + 9*T^2)/(T + 16)^2 ]]>,
<[11,1],53,[[ 28*T^5 + T^4 + 23*T^3 + 40*T^2 + 15*T ], [ (49*T^6 +
44*T^5 + 38*T^4 )/(T^2 + 42*T + 5)^2 ]]>
];

// The formula for Delta_p in Lemma 4.1

function Delta_p(f)
  ff := Factorization(f);
  g := &*[m[1]^m[2] : m in ff | IsCyclotomicPolynomial(m[1])];
  h := ExactQuotient(f,g);
  D := Evaluate(h,1)*Evaluate(h,-1);
  return PowerFreePart(D,2);
end function;

// The number of irreducible components of multiplicity one
// (as appears in Lemma 4.2)

function GeometricTamagawaNumber(kod)
  kod := Sprint(kod);
  if kod in ["I0","II","II*"] then return 1; end if;
  if kod in ["III","III*"] then return 2; end if;
  if kod in ["IV","IV*"] then return 3; end if;
  if kod[#kod] ne "*" then
    n := StringToInteger(&cat[kod[i] : i in [2..#kod]]);
    assert kod eq Sprintf("I%o",n);
    return n;
  else
    n := StringToInteger(&cat[kod[i] : i in [2..#kod-1]]);
    assert kod eq Sprintf("I%o*",n);
    return 4;
  end if;
end function;

TamagawaProduct := func<E|
    &*[GeometricTamagawaNumber(x[5])^Degree(x[1])
                                : x in LocalInformation(E)]>;
BadFibres := func<E|
  {* Sprint(x[5])^^Degree(x[1]) : x in LocalInformation(E) *} >;

TopologicalEulerNumber := func<E|
    &+[Degree(x[1])*x[2] : x in LocalInformation(E) ] >;

print "Checking the polynomials f_p in Table 5";

maxr := 2;
// should be 3 for a complete verification (but then takes several minutes)
print "We compute the n_r for r <=",maxr;

PP<X> := PowerSeriesRing(Rationals());
R2<u,v> := PolynomialRing(Rationals(),2);
P3<[z]> := ProjectiveSpace(Rationals(),3);

for dd in Table5 do
  N,eps := Explode(dd[1]);
  E := Z(N,eps);
  printf "(N,eps) = (%o,%o)\n",N,eps;
  assert TopologicalEulerNumber(E) in [24,36];  
  assert exists(xx){x[2] : x in Table3 | x[1] eq [N,eps]};
  bqfs := [];
  for ee in dd[2] do
    p,fp,D := Explode(ee);
    printf "p = %o  ",p;
    assert Degree(fp) in [22,34]; // deg(f_p) = b_2
    assert Evaluate(fp,0) eq 1;
    assert x^Degree(fp)*Evaluate(fp,1/x) in {fp,-fp}; // functional eqn
    k<T> := FunctionField(GF(p));
    Ep := EllipticCurve([k!a : a in Eltseq(E)]);
    xxp := xx;
    if [N,eps] notin [[9,1],[12,1]] then 
      for x in Table4 do
        if x[1] eq [N,eps] and LegendreSymbol(x[2],p) eq 1 then
          Append(~xxp,x[3]);
        end if;
      end for;
    end if;
    MWgens := [Points(Ep,x)[1]: x in xxp];
    assert Regulator(MWgens) ne 0;
    Frobp := FrobeniusActionOnTrivialLattice(Ep);
    cp := P!CharacteristicPolynomial(Frobp);
    cp *:= (x - 1)^#MWgens;
    fp1 := ExactQuotient(fp,cp);
    assert Degree(fp1) le 6;
    // So we only need to compute n_r for r <= 3
    // (N.B. in our examples with Degree(fp1) = 6
    // the middle coefficient is non-zero, so we know
    // the sign in the functional equation is +)
    nn := [];
    for r in [1..maxr] do
      nr := NumberOfPointsOnSurface(Ep,r);
      printf "n%o = %o  ",r,nr;
      Append(~nn,nr);
    end for;
    print "";
    zeta := Exp(&+[(nn[r] - p^(2*r) - 1)/r*X^r : r in [1..maxr]]);
    Fp := Evaluate(1/zeta,X/p) + O(X^(maxr+1));
    assert Fp eq PP!fp + O(X^(maxr+1));
    assert Delta_p(fp) eq D; // a check using Lemma 4.1
    if [N,eps] in [[9,2],[10,1],[10,3],[11,1]] then
      assert exists(QQ){x[3] : x in ExtraGenerators
                     | x[1] eq [N,eps] and x[2] eq p };
      QQ := [#Q eq 1 select Points(Ep,Q[1])[1] else Ep!Q : Q in QQ];
      DD := TamagawaProduct(Ep)*Regulator(MWgens cat QQ);
      assert PowerFreePart(DD,2) eq D; // a check using Lemma 4.2
      rhop := #(MWgens cat QQ) + 2;
      for x in LocalInformation(Ep) do
        pl,val,f,c,kod := Explode(x);
        m := val + 1 - f; // Ogg's formula
        rhop +:= Degree(pl)*(m - 1); // Shioda-Tate formula
      end for;
      assert rhop eq 30;
    end if;
    if [N,eps] in [[10,1],[10,3],[11,1]] then
      assert BadFibres(E) eq BadFibres(Ep);
      a := Regulator(MWgens cat [QQ[1]]);
      bb := Regulator(MWgens cat [QQ[1]+QQ[2]]);
      c := Regulator(MWgens cat [QQ[2]]);
      Append(~bqfs,a*u^2 + (bb - a - c)*u*v + c*v^2);
    end if;  
  end for;
  if #bqfs eq 2 then
    cc := [RationalGCD(Coefficients(bqf)): bqf in bqfs];
    printf "Quadratic form : (%o)*(%o) = ",cc[1],bqfs[1]/cc[1];
    printf "(%o)*(%o)\n",cc[2],bqfs[2]/cc[2];
    LHS := Evaluate(bqfs[1],[z[1],z[2]]);
    RHS := Evaluate(bqfs[2],[z[3],z[4]]);   
    V := Scheme(P3,LHS - RHS);
    printf "Not locally soluble at p = ";
    print [p : p in [2..30] |
               IsPrime(p) and not IsLocallySoluble(V,p)];
  end if;
end for;       

////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////

print "";
print "Some examples of pairs of N-congruent elliptic curves";
print "";

// Finally we use the formulae in the paper to construct some
// pairs of N-congruent elliptic curves over Q. Such curves have
// traces of Frobenius (at all good primes) that are congruent mod N.
// This provides some check on our calculations.

// The following are some test points on the surfaces Z(N,eps).
// Each surface is defined by a Weierstrass equation
// y^2 + a1(T)*x*y + a3(T)*y = x^3 + a2(T)*x^2 + a4(T)*x + a6(T)
// and we specify points as [T,x,y].

// These points were chosen at random, avoiding the finitely
// many curves on each surface where our formulae do not apply.

test_points := [
<[6,5], [ [ -1, 376, -8836 ], [ 3/2, 33/8, -121/16 ],
      [ -4, 5005, -511225 ], [ 4/3, 80/81, 320/729 ],
      [ -3/2, 129/8, -1849/16 ] ]>,
<[7,3], [ [ -3/2, 144, -336 ], [ -2, 208, 468 ],
  [ -13/6, 4/9, -3068/81 ], [ 1, 49, -21 ], [ 7/2, 144, 3024 ] ]>,
<[8,3], [ [ 2, 5, 132 ], [ 1/3, 427/27, 2057/27 ], [ -2, 5, -132 ],
   [ -3/2, 378, 29403/4 ], [ -2/3, 4316/27, -55556/27 ] ]>,
<[8,5], [ [ 3/2, 5/2, 75/2 ], [ 11/4, 1225/32, 4375/64 ],
   [ 4, 735, 18900 ], [ 0, 63, 420 ], [ 3, 16, -24 ] ]>,
<[8,7], [ [ -3/2, 35/8, -315/16 ], [ -1/2, 81/64, -1215/512 ],
   [ -2, 69, -1380 ], [ 2, 9, -360 ], [ -3/4, -119/1024, 5355/16384 ] ]>,
<[9,1], [ [ -3/4, 9/64, -63/1024 ], [ -3/2, 63/4, -441/4 ],
  [ -5/2, 475/4, -9025/4 ], [ 1, 22, -44 ], [ 3/2, -9/2, -225/4 ] ]>,
<[9,2], [ [ 1/3, -448/243, 2048/729 ], [ 7, 31104, -16702848 ],
    [ 3/2, 35/16, -5501/128 ], [ -5/4, 2187/512, 7814151/16384 ],
    [ -5, -5184, 2985984 ] ]>,
<[10,1], [ [ 3, -105444/25, 32265864/125 ], [ -1, 850, 12750 ],
    [ 4/3, 640/81, 102400/2187 ], [ 4, -40320, 2764800 ],
    [ 1/2, 1185/1024, -17775/32768 ] ]>,
<[10,3], [ [ 3, 32, -5248 ], [ 4, 48708/49, 11560032/343 ],
    [ -5/2, 23213/8, 2065957/8 ], [ -1/2, 225/8, 7875/32 ],
    [ -4/3, -1496/81, -32912/729 ] ]>,   
<[11,1], [ [ 3, 121, -363 ], [ 5/6, 245/144, -5845/31104 ],
  [ -1, 15, -75 ], [ 7/2, 621/16, -1863/128 ], [ -3, 1325, 111035 ] ] >,
<[12,1], [ [ -5, 0, 0 ], [ -3, -1280, 12800 ],
  [ 1/3, -560/81, -172480/729 ], [ 3/2, -455/4, -15925/16 ],
  [ 7, 26000, 1352000 ] ] > ];

for Neps in [[6,5],[7,3],[8,3],[8,5],[8,7],[9,1],[9,2],
             [10,1],[10,3],[11,1],[12,1]] do
  N,eps := Explode(Neps);
  printf "Some examples in the case (N,eps) = (%o,%o)\n",N,eps;  
  trfrob := func<E|Vector(Integers(N),TracesOfFrobenius(E,1000))>;
  assert exists(pts){ t[2] : t in test_points | t[1] eq [N,eps] };
  if N notin [6,10] then 
    E,aa,xx := ModularInterpretation(N,eps);
  else
    E,uu,yy := MapToZstar(N,eps);
  end if;
  for pt in pts do
    T0,x0,y0 := Explode(pt);
    ET := EllipticCurve([Evaluate(a,T0): a in Eltseq(E)]);
    printf "T = %o,  P = %o\n",T0,ET![x0,y0];
    if N notin [6,10] then 
      a,b := Explode([MyEvaluate(f,pt): f in aa]);
      xx0 := [MyEvaluate(f,pt): f in xx];
      E1 := EllipticCurve([a,b]);
      E2 := CongruentEllipticCurve(N,eps,[a,b],xx0);
    else
      u,v,w,y := Explode([MyEvaluate(f,pt): f in uu cat [yy]]);
      E1,E2 := CongruentEllipticCurves(N,eps,[u,v,w],y);
    end if;
    if N eq 11 then
    // In this case there are two candidates for E2 (quadratic twists by -1)
      EE2 := E2;
      E2 := EE2[1];
      diffr := trfrob(E1) - trfrob(E2);
      if #[x : x in Eltseq(diffr) | x ne 0] gt 5 then
        E2 := EE2[2];
      end if;
    end if;
    printf "N(E1) = %o  N(E2) = %o\n",Conductor(E1),Conductor(E2);
    diffr := trfrob(E1) - trfrob(E2);
    assert #[x : x in Eltseq(diffr) | x ne 0] lt 5;
  end for;
  print "";
end for;

