
// 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-formulae.m" we record some of the formulae
// relating to the the modular diagonal quotient surfaces Z(N,eps),
// and the twisted modular curves X_E(N,eps). These formulae are
// needed to run the other file "congr-checkpaper.m".

Deriv := Derivative;
MC := MonomialCoefficient;
MD := MonomialsOfDegree;
MyEvaluate := func<f,pt|
  eval Sprintf("%o where T,x,y is Explode(%o)",f,pt)>;
Prim := func<y|[x/d: x in xx] where d is RationalGCD(xx) 
                              where xx is Eltseq(y)>;
K<T> := FunctionField(Rationals());

// The elliptic surfaces in Theorems 1.1 and 1.2.

function ModularDiagonalQuotientSurface(N,eps)
  case [N,eps] : 
  when [6,5] :
    aa := [3*T*(T - 2),-6*(T - 1)*(T^3 - 2),
           2*(T - 1)*(T + 2)^2*(T^3 - 2),0,0];
  when [7,3] : 
    aa := [0,4*T^4 + 4*T^3 - 51*T^2 - 2*T - 50,0,
           (6*T + 25)*(52*T^2 - 4*T + 25),0];
  when [8,3] : 	   
    aa := [0,-(3*T^2 - 7),0,-4*T^2*(4*T^4 - 15),
           4*T^2*(53*T^4 + 81*T^2 + 162)];
  when [8,5] :	   
    aa := [0,-2*(T^2 + 19),0,
           -(4*T^2 - 49)*(T^4 - 6*T^2 + 25),0];
  when [9,1] :
    aa := [6*T^2 + 3*T + 2,-(16*T^4 + 12*T^3 + 9*T^2 + 6*T + 1),
           T^2*(T + 1)*(4*T^3 + 9*T + 9),0,0];
  when [12,1] : 	   
    aa := [2*(5*T^2 + 9),(T^2 + 3)*(11*T^2 + 1),
           96*(T^2 + 3)*(T^2 + 1)^2,0,0];
  when [8,7] : 	   
    aa := [0,2*(4*T^6 - 15*T^4 + 14*T^2 - 1),0,
           (T^2 - 1)^4*(16*T^4 - 24*T^2 + 1),0];
  when [9,2] : 	   
    aa := [3*(4*T^3 + T^2 - 2),- 3*(T + 1)*(T^3 - 1)*(9*T^2 + 2*T + 1),
           (T - 1)^3*(T^3 - 1)*(4*T^3 - 3*T - 7),0,0];
  when [10,1] :
    aa := [-(3*T - 2)*(6*T^2 - 5*T - 2),T^2*(T - 1)
           *(27*T^3 - 54*T^2 + 16*T + 12), -4*T^2*(T - 1)^2
	   *(4*T^2 - 2*T - 1)*(27*T^3 - 54*T^2 + 16*T + 12),0,0];
  when [10,3] : 
    aa := [T^3 - 8*T^2 - 9*T - 8,2*(3*T + 2)*(T^3 - T^2 - 3*T - 3),
           2*T^2*(T^3 - T^2 - 3*T - 3)*(7*T^2 + 2*T + 3),0,0];
  when [11,1] : 	   
    aa := [T^3 + T,-(4*T^5 - 17*T^4 + 30*T^3 - 18*T^2 + 4),
           0,T^2*(2*T - 1)*(3*T^2 - 7*T + 5)^2,0];
  end case;
  return EllipticCurve(aa);
end function;

// Equations for the twisted modular curves X_E(N,eps).
// (Taken from earlier papers, as referenced below)

function XE_Equations(N,eps,aa)
  R := Universe(aa);
  if R cmpeq Integers() then R := Rationals(); end if;
  a,b := Explode(aa);
  case [N,eps] : 
  when [7,3] :
    RR<x,y,z> := PolynomialRing(R,3);
    // Equation for X_E(7,3) --- Poonen, Schaefer, Stoll.
    FF := [ -a^2*x^4 + 2*a*b*x^3*y - 12*b*x^3*z - 6*(a^3 + 6*b^2)*x^2*y^2 
       + 6*a*x^2*z^2 + 2*a^2*b*x*y^3 - 12*a*b*x*y^2*z + 18*b*x*y*z^2 
       + (3*a^4 + 19*a*b^2)*y^4 - 2*(4*a^3 + 21*b^2)*y^3*z 
       + 6*a^2*y^2*z^2 - 8*a*y*z^3 + 3*z^4 ];
  when [8,3] :
    RR<x0,x1,x2,x3,x4> := PolynomialRing(R,5);
    // Equations for X_E(8,3) --- Chen
    FF := [ x0^2 + 2*b*x1*x3 + a*x2^2 - 6*b*x2*x4 - a^2*x4^2,
       2*x0*x1 + 2*a*x1*x3 + 4*a*x2*x4 - b*x3^2 - 18*b*x4^2,
       2*x0*x3 - x1^2 + x2^2 + a*x3^2 + 3*a*x4^2 ];
  when [8,5] : 
    RR<x0,x1,x2,x3,x4> := PolynomialRing(R,5);
    // Equations for X_E(8,5) --- Chen
    FF := [ -a*x3^2 + 2*x1*x3 + x2^2 - (5*a*x4^2 - 3*x0^2),
       -2*a*x2*x3 - b*x3^2 + 2*x1*x2 - (2*a*x0*x4 - 6*b*x4^2),
       -2*b*x2*x3 + x1^2 - (4*a^2*x4^2 - 4*a*x0^2 - 6*b*x0*x4) ];
  when [8,7] :
    RR<x0,x1,x2,x3,x4> := PolynomialRing(R,5);
    // Equations for X_E(8,7) --- Chen
    FF := [ 3*x0^2 + a*x4^2 - a*x3^2 + 2*x1*x3 + x2^2,
       4*a*x0*x4 + 6*b*x4^2 - 2*a*x2*x3 - b*x3^2 + 2*x1*x2,
       a*x0^2 + 6*b*x0*x4 - a^2*x4^2 - 2*b*x2*x3 + x1^2 ];
  when [9,1] :
    RR<x,y,z,t> := PolynomialRing(R,4);
    // Equations for X_E(9,1) --- Fisher
    FF := [ x^2*t + 6*x*y*z + 6*b*x*t^2 + 6*y^3 - 9*a*y^2*t + 6*a^2*y*t^2 
       - 3*b*z^3 + 3*a^2*z^2*t + 9*a*b*z*t^2  - (a^3 - 12*b^2)*t^3,
       x^2*z + 6*x*y^2 - 6*a*x*y*t + 2*a^2*x*t^2 - 9*a*y^2*z - 18*b*y*z^2 
       + 12*a^2*y*z*t + a^2*z^3 + 9*a*b*z^2*t - 3*a^3*z*t^2 + a^2*b*t^3 ];
  when [9,2] :
    RR<x,y,z,t> := PolynomialRing(R,4);
    // Equations for X_E(9,2) --- Fisher
    FF := [ 9*x^2*y + 3*x^2*z - 6*a*x*y*t + 6*b*x*t^2 - 6*a*y^3 + 27*b*y^2*t 
       + 3*a*y*z^2 + 18*b*y*z*t + 3*a^2*y*t^2 + a*z^3 + 3*b*z^2*t 
       + a^2*z*t^2 - a*b*t^3,
       x^3 + 6*a*x*y*z + 18*b*x*y*t + 3*a*x*z^2 + 6*b*x*z*t + a^2*x*t^2 
       + 9*b*y^3 + 6*a^2*y^2*t - 9*b*y*z^2 + 6*a^2*y*z*t - 3*a*b*y*t^2 
       - 4*b*z^3 + 2*a^2*z^2*t + 2*b^2*t^3 ];
  when [11,1] :
    RR<v,w,x,y,z> := PolynomialRing(R,5);
    // Cubic form determining X_E(11,1) --- Fisher
    FF := [ v^3 + a*v^2*z - 2*a*v*x^2 + 4*a*v*x*y - 6*b*v*x*z + a*v*y^2 
       + 6*b*v*y*z + a^2*v*z^2 - w^3 + a*w^2*z - 4*a*w*x^2 - 12*b*w*x*z 
       + a^2*w*z^2 - 2*b*x^3 + 3*b*x^2*y + 2*a^2*x^2*z + 6*b*x*y^2 
       + 4*a*b*x*z^2 + b*y^3 - a^2*y^2*z + a*b*y*z^2 + 2*b^2*z^3 ];
  when [12,1] :      
    RR<x,y,z,u0,u1,u2> := PolynomialRing(R,6);
    // Equations for X_E(12,1) --- Chen
    FF := [ -x^2*z + a*x*y^2 + 6*b*y^3 - 6*a*y^2*z - 12*z^3,
       x^2 + 12*x*z + 36*z^2 - 2*u0*u2 - u1^2 + a*u2^2,
       4*a*x*y + 36*b*y^2 - 24*a*y*z - 2*u0*u1 + 2*a*u1*u2 + b*u2^2,
       8*a*x*z - 4*a^2*y^2 - u0^2 + 2*b*u1*u2 ];
  end case;
  return FF;
end function;

// Some formulae in the case N = 9 taken from [10]

function QuarticForm(c4,c6,FF,sign)
  P<T> := PolynomialRing(Universe(FF));
  HH := [Matrix(4,4,[Deriv(Deriv(F,i),j): i,j in [1..4]]): F in FF];
  DD := Determinant(HH[1] + T*HH[2]);
  if sign eq 1 then
    D := (1/144)*Coefficient(DD,0);
    qq := 144*(1 - 54*c4*T^2 + 216*c6*T^3 - 243*c4^2*T^4);
  else
    if c4 ne 0 then
      D := (1/c4)*Coefficient(DD,0);
    else
      D := -1/(12*c6)*Coefficient(DD,1);
    end if;
    qq := c4 - 12*c6*T + 54*c4^2*T^2
          - 108*c4*c6*T^3 - (243*c4^3-324*c6^2)*T^4;
  end if;
  assert DD eq qq*D;
  return D;
end function;

function MapFromX9TwistToP1(F,pt)
  P<X> := PolynomialRing(Rationals());
  Jmat := Matrix(4,2,[Evaluate(Deriv(f,i),pt): f in F,i in [1..4]]);
  km := KernelMatrix(Jmat);
  km := LLL(Matrix(Basis(PureLattice(Lattice(km)))));
  mat := Matrix(2,4,[pt,Eltseq(km[1])]);
  if Rank(mat) lt 2 then 
    mat := Matrix(2,4,[pt,Eltseq(km[2])]);
  end if;
  subst := [mat[1,i] + X*mat[2,i]: i in [1..4]];
  aa := [Coefficient(Evaluate(f,subst),2): f in F];
  return Prim([aa[2],-aa[1]]);
end function;

// Some formulae in the case N = 11 taken from [9]

function c4Invariant(F)
  R := Parent(F);
  PP<T> := PolynomialRing(R);
  KK<t> := quo<PP|T^4>;
  Hmat1 := Matrix(5,5,[Deriv(Deriv(F,i),j):i,j in [1..5]]);
  H := (1/32)*Determinant(Hmat1); 
  Hmat2 := Matrix(5,5,[Deriv(Deriv(H,i),j): i,j in [1..5]]);
  poly := Determinant(Hmat1 + t*Hmat2);
  assert Coefficient(poly,0) eq 32*H;
  return (1/8)*Coefficient(poly,3);
end function;

function QuinticForm(inv)
  c4,DD := Explode(inv);
  R<v,w,x,y,z> := PolynomialRing(Universe(inv),5);
  Q := 78096*DD*v^4*w + 743*DD*c4^2*v^4*x + 4816*DD*c4*v^4*y +
    3254*DD*c4^2*v^4*z + 84660*DD*v^3*w^2 + 919*DD*c4^2*v^3*w*x -
    58228*DD*c4*v^3*w*y - 5961*DD*c4^2*v^3*w*z - 31725*DD^2*c4*v^3*x^2
    - 402192*DD^2*v^3*x*y + 147816*DD^2*c4*v^3*x*z +
    6040*DD*c4^2*v^3*y^2 + (3317472*DD^2 - 1309*DD*c4^3)*v^3*y*z +
    (206244*DD^2*c4 - 276*DD*c4^4)*v^3*z^2 + 848748*DD*v^2*w^3 +
    10749*DD*c4^2*v^2*w^2*x - 208728*DD*c4*v^2*w^2*y +
    95511*DD*c4^2*v^2*w^2*z - 116055*DD^2*c4*v^2*w*x^2 + (-173232*DD^2
    - 919*DD*c4^3)*v^2*w*x*y - 1349352*DD^2*c4*v^2*w*x*z +
    39520*DD*c4^2*v^2*w*y^2 + (3297888*DD^2 - 8207*DD*c4^3)*v^2*w*y*z
    + (-8068356*DD^2*c4 - 317*DD*c4^4)*v^2*w*z^2 +
    1095648*DD^3*v^2*x^3 - 10447*DD^2*c4^2*v^2*x^2*y + (15372288*DD^3
    + 8026*DD^2*c4^3)*v^2*x^2*z - 268080*DD^2*c4*v^2*x*y^2 +
    109224*DD^2*c4^2*v^2*x*y*z + (132860736*DD^3 +
    39537*DD^2*c4^3)*v^2*x*z^2 + (-3187968*DD^2 -
    3806*DD*c4^3)*v^2*y^3 + (880128*DD^2*c4 + 1401*DD*c4^4)*v^2*y^2*z
    + 15444*DD^2*c4^2*v^2*y*z^2 + (395788032*DD^3 +
    96126*DD^2*c4^3)*v^2*z^3 - 1200228*DD*v*w^4 + 3924*DD*c4^2*v*w^3*x
    + 611412*DD*c4*v*w^3*y - 270767*DD*c4^2*v*w^3*z -
    236205*DD^2*c4*v*w^2*x^2 + (-8028720*DD^2 -
    3630*DD*c4^3)*v*w^2*x*y + 700272*DD^2*c4*v*w^2*x*z -
    105870*DD*c4^2*v*w^2*y^2 + (-49864032*DD^2 +
    64964*DD*c4^3)*v*w^2*y*z + (9933516*DD^2*c4 -
    15154*DD*c4^4)*v*w^2*z^2 + (135696*DD^3 - 1072*DD^2*c4^3)*v*w*x^3
    + 34582*DD^2*c4^2*v*w*x^2*y + (-683424*DD^3 -
    10408*DD^2*c4^3)*v*w*x^2*z + 886848*DD^2*c4*v*w*x*y^2 -
    264240*DD^2*c4^2*v*w*x*y*z + (-11943936*DD^3 +
    78489*DD^2*c4^3)*v*w*x*z^2 + (-13607040*DD^2 -
    4608*DD*c4^3)*v*w*y^3 + (19470816*DD^2*c4 - 848*DD*c4^4)*v*w*y^2*z
    + (-6935616*DD^2*c4^2 + 317*DD*c4^5)*v*w*y*z^2 + (-497539584*DD^3
    + 924300*DD^2*c4^3)*v*w*z^3 + 4088*DD^3*c4^2*v*x^4 -
    304120*DD^3*c4*v*x^3*y - 8466*DD^3*c4^2*v*x^3*z + (-1218528*DD^3 +
    4490*DD^2*c4^3)*v*x^2*y^2 + (-6074784*DD^3*c4 -
    1204*DD^2*c4^4)*v*x^2*y*z + 297810*DD^3*c4^2*v*x^2*z^2 +
    63640*DD^2*c4^2*v*x*y^3 + (-108570240*DD^3 -
    3538*DD^2*c4^3)*v*x*y^2*z + (22483008*DD^3*c4 +
    2955*DD^2*c4^4)*v*x*y*z^2 + (-7726968*DD^3*c4^2 -
    2248*DD^2*c4^5)*v*x*z^3 + (1936480*DD^2*c4 + 927*DD*c4^4)*v*y^4 +
    (-1324560*DD^2*c4^2 - 467*DD*c4^5)*v*y^3*z + (-305731584*DD^3 +
    323442*DD^2*c4^3)*v*y^2*z^2 + (181723392*DD^3*c4 +
    799*DD^2*c4^4)*v*y*z^3 + (-53461728*DD^3*c4^2 -
    15542*DD^2*c4^5)*v*z^4 + (2662068*DD + 1815*c4^3)*w^5 -
    22829*DD*c4^2*w^4*x + 47352*DD*c4*w^4*y - 48793*DD*c4^2*w^4*z +
    548049*DD^2*c4*w^3*x^2 + (-18046800*DD^2 - 11102*DD*c4^3)*w^3*x*y
    + (8968896*DD^2*c4 + 1222*DD*c4^4)*w^3*x*z - 74634*DD*c4^2*w^3*y^2
    + (-103691232*DD^2 - 77590*DD*c4^3)*w^3*y*z + (37026036*DD^2*c4 +
    23002*DD*c4^4)*w^3*z^2 - 3794112*DD^3*w^2*x^3 -
    76569*DD^2*c4^2*w^2*x^2*y + (-76011264*DD^3 +
    34360*DD^2*c4^3)*w^2*x^2*z + (1160640*DD^2*c4 +
    1815*DD*c4^4)*w^2*x*y^2 + 215352*DD^2*c4^2*w^2*x*y*z +
    (-399940416*DD^3 + 182544*DD^2*c4^3)*w^2*x*z^2 + (-268416*DD^2 +
    27312*DD*c4^3)*w^2*y^3 + (8347968*DD^2*c4 - 87*DD*c4^4)*w^2*y^2*z
    + 5276556*DD^2*c4^2*w^2*y*z^2 + (-1570772736*DD^3 -
    1524066*DD^2*c4^3)*w^2*z^3 - 4704*DD^3*c4^2*w*x^4 +
    481768*DD^3*c4*w*x^3*y - 399014*DD^3*c4^2*w*x^3*z + (28550304*DD^3
    + 16365*DD^2*c4^3)*w*x^2*y^2 + (-12310848*DD^3*c4 -
    14946*DD^2*c4^4)*w*x^2*y*z + (-3905442*DD^3*c4^2 +
    1213*DD^2*c4^5)*w*x^2*z^2 + 167736*DD^2*c4^2*w*x*y^3 +
    (506307456*DD^3 + 91410*DD^2*c4^3)*w*x*y^2*z + (-304356096*DD^3*c4
    - 161654*DD^2*c4^4)*w*x*y*z^2 + (3995784*DD^3*c4^2 +
    12626*DD^2*c4^5)*w*x*z^3 + 5196896*DD^2*c4*w*y^4 -
    3720112*DD^2*c4^2*w*y^3*z + (2376511488*DD^3 +
    1293222*DD^2*c4^3)*w*y^2*z^2 + (-1450455552*DD^3*c4 -
    588877*DD^2*c4^4)*w*y*z^3 + (130012992*DD^3*c4^2 +
    49798*DD^2*c4^5)*w*z^4 - 1920*DD^4*c4*x^5 + (1790208*DD^4 +
    1272*DD^3*c4^3)*x^4*y + 633024*DD^4*c4*x^4*z -
    3476*DD^3*c4^2*x^3*y^2 + (140363136*DD^4 +
    113974*DD^3*c4^3)*x^3*y*z + (6337296*DD^4*c4 -
    4469*DD^3*c4^4)*x^3*z^2 + (-4372512*DD^3*c4 -
    2762*DD^2*c4^4)*x^2*y^3 + (1425624*DD^3*c4^2 +
    602*DD^2*c4^5)*x^2*y^2*z + (2289731328*DD^4 +
    1689870*DD^3*c4^3)*x^2*y*z^2 + (-28416096*DD^4*c4 -
    103935*DD^3*c4^4)*x^2*z^3 + (-47748096*DD^3 -
    40052*DD^2*c4^3)*x*y^4 + (-49195008*DD^3*c4 -
    24132*DD^2*c4^4)*x*y^3*z + (29284560*DD^3*c4^2 +
    11512*DD^2*c4^5)*x*y^2*z^2 + (19383017472*DD^4 +
    10628208*DD^3*c4^3)*x*y*z^3 + (-1153336320*DD^4*c4 -
    813060*DD^3*c4^4)*x*z^4 - 244000*DD^2*c4^2*y^5 + (-458422272*DD^3
    - 62696*DD^2*c4^3)*y^4*z + (-92943360*DD^3*c4 -
    120506*DD^2*c4^4)*y^3*z^2 + (108989280*DD^3*c4^2 +
    54983*DD^2*c4^5)*y^2*z^3 + (59808264192*DD^4 + 25465680*DD^3*c4^3
    - 1596*DD^2*c4^6)*y*z^4 + (-6097006080*DD^4*c4 -
    2940300*DD^3*c4^4)*z^5;
  return Q;
end function;

// Given a point xx on X_E(N,eps) where E is the elliptic curve
// y^2 = x^3 + aa[1]*x + aa[2], we compute the corresponding
// elliptic curve N-congruent to E with power eps.
// References are the same as for "XE_Equations" as above,
// except in the case (N,eps) = (7,3), for which see [9].

function CongruentEllipticCurve(N,eps,aa,xx)
  FF := XE_Equations(N,eps,aa);
  if N eq 11 then
    F := FF[1];
    Hmat := Matrix(5,5,[Deriv(Deriv(F,i),j):i,j in [1..5]]);
    FF := Minors(Hmat,4);
  end if;
  assert forall{ F : F in FF | Evaluate(F,xx) eq 0 };
  a,b := Explode(aa);
  Delta := -(4*a^3 + 27*b^2);
  case [N,eps] :
  when [7,3] :
    Hmat := Matrix(3,3,[Deriv(Deriv(FF[1],i),j): i,j in [1..3]]);
    Psi_6 := (-1/54)*Determinant(Hmat);
    vv := Matrix(3,1,[Deriv(Psi_6,i): i in [1..3]]);
    HHmat := HorizontalJoin(Hmat,vv);
    HHmat := VerticalJoin(HHmat,Matrix(1,4,Eltseq(vv) cat [0]));
    Psi_14 := (1/9)*Determinant(HHmat);
    Jmat := Matrix(3,3,[Deriv(f,i): i in [1..3],f in [FF[1],Psi_6,Psi_14]]);
    Psi_21 := (1/14)*Determinant(Jmat);
    x,y,z := Explode(xx);
    dd := -Delta*(7*a*x^2*y - 6*x^2*z - 3*a^2*y^3 + 8*a*y^2*z - 3*y*z^2);
    cc4 := Evaluate(Psi_14,xx);
    cc6 := Evaluate(Psi_21,xx);
    E := EllipticCurve([-27*dd^2*cc4,-54*dd^3*cc6]);
  when [8,3] : 
    x0,x1,x2,x3,x4 := Explode(xx);
    t := (-3*b*x2 - 2*a^2*x4)/(6*a*x2 - 27*b*x4);
    _,cc4,cc6 := HessePolynomials(4,1,[-a/27,-b/54]:Variables := [t,1]);
    E := EllipticCurve([-27*Delta^2*cc4,-54*Delta^3*cc6]);
  when [8,5] : 
    x0,x1,x2,x3,x4 := Explode(xx);
    t := x0/(3*x4);  
    _,cc4,cc6 := HessePolynomials(4,1,[-a/27,-b/54]:Variables := [t,1]);
    E := EllipticCurve([-27*cc4,-54*cc6]);
  when [8,7] : 
    x0,x1,x2,x3,x4 := Explode(xx);
    t := x0/(3*x4);  
    _,cc4,cc6 := HessePolynomials(4,1,[-a/27,-b/54]:Variables := [t,1]);
    E := EllipticCurve([-27*Delta^2*cc4,-54*Delta^3*cc6]);
  when [9,1] : 
    D := QuarticForm(-a/27,-b/54,FF,1);
    assert Evaluate(D,xx) ne 0;  
    aa := MapFromX9TwistToP1(FF,xx);
    aa := Prim([aa[1],-3*aa[2]]);
    _,cc4,cc6 := HessePolynomials(3,1,[-a/27,-b/54]:Variables := aa);
    E := EllipticCurve([-27*cc4,-54*cc6]);
  when [9,2] : 
    D := QuarticForm(-a/27,-b/54,FF,-1);
    assert Evaluate(D,xx) ne 0;
    c4,c6 := Explode([-a/27,-b/54]);
    aa := MapFromX9TwistToP1(FF,xx);
    aa := Prim([c6*aa[1] - 3*c4^2*aa[2],-c4*aa[1] + 3*c6*aa[2]]);
    _,cc4,cc6 := HessePolynomials(3,-1,[c4,c6]:Variables := aa);
    E := EllipticCurve([-27*2^2*cc4,-54*2^3*cc6]);
  when [11,1] :
    v,w,x,y,z := Explode(xx);
    d := Evaluate(QuinticForm([-a/27,-16*Delta/6^12]),
      [ -9*b*v - 18*b*w - 4*a^2*x + 2*a^2*y - 3*a*b*z, 9*b*v - 4*a^2*x
       - 2*a^2*y - 3*a*b*z, 419904*x, 648*a*x + 972*b*z, 34992*y ]);
    assert b*d ne 0; 
    mu := PowerFreePart(d/(2*3^3*5*11^2*b),4);
    cc4 := mu*Evaluate(c4Invariant(F),xx);
    DD := -2^16*mu^3*Delta^16*Evaluate(F,xx)^11;
    flag,cc6 := IsSquare(cc4^3 - 1728*DD); 
    assert flag;
    E := [EllipticCurve([-27*cc4,54*k*cc6]): k in [-1,1]];
  when [12,1] : 
    X,Y,Z,u0,u1,u2 := Explode([xx[i]/xx[3] : i in [1..6]]);
    Delta *:= 16;
    x := a*X + 6*b*Y - 2*a;
    y := a^2*X*Y + 6*b*X + 6*a*b*Y^2 - 6*a^2*Y;
    la := (x^3*y - 108*b*x^3 - 8*Delta*y)/(18*(x^4 + 12*a*x^3 + 4*Delta*x));
    _,cc4,cc6 := HessePolynomials(3,1,[-a/27,-b/54]:Variables := [la,1]);
    E := EllipticCurve([-27*cc4,-54*cc6]);
  end case;
  return E;
end function;

// Each point on Z(N,eps) determines a pair of elliptic curves
// E1 and E2. We compute tuples of rational functions aa and xx
// such that E1 has equation y^2 = x^3 + aa[1]*x + aa[2], and
// xx is the point corresponding to E2 on X_E1(N,eps).

// The formulae used are based on those in Section 2 (or [6],[10]
// in the cases N = 8,9)

function ModularInterpretation(N,eps)
  E := ModularDiagonalQuotientSurface(N,eps);
  KE<x,y> := FunctionField(E);
  case [N,eps] :
  when [7,3] : 
    R<x1,x2> := PolynomialRing(K,2);
    T0 := (6*T - 3)/(4*T + 4);
    c := (T0^2 + 1)*(3*T0^2 - 8*T0 + 3);
    h := T0^3*(19*T0 - 42)*x1^2 + 2*T0*(T0 - 3)^2*x1*x2 - 6*(T0^2 - 1)*x2^2;
    f := 36*T0^2*x1^2*x2^2 - 2*(T0 - 6)*x1*x2^3 + x2^4;
    R<x> := PolynomialRing(K);
    C := HyperellipticCurve(Evaluate(h^2 + 4*c*f,[1,x]));
    E1,map1 := EllipticCurve(C,C![0,T0^3*(19*T0 - 42)]);
    flag,map2 := IsIsomorphic(E,E1);
    assert flag;
    mapEtoC := Expand(map2 * map1^(-1));
    ff := DefiningPolynomials(mapEtoC);
    x2 := KE!(ff[1]/ff[3]);
    y := KE!(ff[2]/ff[3]^2);
    y := (y - Evaluate(h,[1,x2]))/(2*c);
    aa := [y,y];
    xx := [x2,T0,y];
  when [8,3] :
    ff := [ (-24*T^11 - 48*T^9 - 198*T^7 + 42*T^5 + 174*T^3 + 54*T)*y
      + (-36*T^6 + 18*T^4 + 18*T^2)*x^2 + (-144*T^10 - 288*T^8 -
      396*T^6 + 504*T^4 + 324*T^2)*x - 144*T^14 - 648*T^12 - 1764*T^10
      - 1782*T^8 - 198*T^6 + 3078*T^4 + 1458*T^2, (-24*T^11 - 48*T^9 -
      198*T^7 + 42*T^5 + 174*T^3 + 54*T)*y + (-72*T^12 - 216*T^10 -
      738*T^8 - 468*T^6 + 648*T^4 + 684*T^2 + 162)*x - 144*T^14 +
      144*T^12 - 324*T^10 + 3816*T^8 + 288*T^6 - 2808*T^4 - 972*T^2,
      (3*T^3*x - 8*T^11 - 12*T^9 - 54*T^7 + 59*T^5 + 63*T^3)*y +
      (-16*T^6 + 25*T^4 + 9*T^2)*x^2 + (8*T^12 - 44*T^10 - 6*T^8 -
      31*T^6 + 280*T^4 + 45*T^2)*x - 48*T^14 - 280*T^12 - 572*T^10 -
      930*T^8 + 228*T^6 + 1026*T^4, (36*T^7 + 18*T^5)*x^2 + (72*T^11 +
      216*T^9 + 18*T^7 - 360*T^5 - 162*T^3)*x - 792*T^11 - 3024*T^9 -
      6498*T^7 - 5508*T^5 - 1458*T^3 ];
    T0 := 1/T;
    v,w,y,z := Explode(ff);
    u := 6*w - 6*v + (6*T0^3 - 6*T0)*z;
    aa := [2*u*v - 3*w^2, u^2*y - 2*u*v*w - T0*u*w*z + 2*w^3 ];
    xx := [-3*T0*u*v - (T0^2 - 1)*u*z + 6*T0*w^2,-3*T0*w,u - 3*w,3*T0,1];
  when [8,5] : 
    ff := [ ((2*T^7 - 16*T^6 + 22*T^5 + 124*T^4 - 388*T^3 + 12*T^2 +
      1104*T - 1080)*x - 4*T^10 + 30*T^9 + 8*T^8 - 632*T^7 + 1568*T^6
      + 2214*T^5 - 14224*T^4 + 13468*T^3 + 26052*T^2 - 62880*T +
      37800)*y + (-T^5 + 10*T^4 - 32*T^3 + 6*T^2 + 156*T - 216)*x^3 +
      (-T^9 + 6*T^8 + 8*T^7 - 106*T^6 + 261*T^5 - 746*T^4 + 1544*T^3 +
      2406*T^2 - 13692*T + 14040)*x^2 + (4*T^12 - 18*T^11 - 112*T^10 +
      605*T^9 + 914*T^8 - 7888*T^7 + 2630*T^6 + 34973*T^5 - 12974*T^4
      - 125152*T^3 + 25338*T^2 + 416580*T - 415800)*x - 4*T^15 +
      20*T^14 + 119*T^13 - 714*T^12 - 1444*T^11 + 11986*T^10 +
      6118*T^9 - 122724*T^8 + 75400*T^7 + 676496*T^6 - 942049*T^5 -
      1942314*T^4 + 4534360*T^3 + 1322250*T^2 - 9586500*T + 6615000,
      ((2*T^4 - 8*T^3 - 2*T^2 + 48*T - 60)*x - 2*T^8 + 8*T^7 + 14*T^6
      - 96*T^5 - 2*T^4 + 488*T^3 - 310*T^2 - 1200*T + 1500)*y +
      (-2*T^6 + 4*T^5 + 20*T^4 - 52*T^3 - 38*T^2 + 168*T - 60)*x^2 +
      (2*T^10 - 38*T^8 + 8*T^7 + 330*T^6 - 96*T^5 - 1634*T^4 + 488*T^3
      + 4540*T^2 - 1200*T - 6000)*x - 4*T^13 + 6*T^12 + 92*T^11 -
      158*T^10 - 992*T^9 + 1848*T^8 + 6424*T^7 - 12796*T^6 - 26020*T^5
      + 55350*T^4 + 63500*T^3 - 143750*T^2 - 75000*T + 187500, (-4*T^4
      + 16*T^3 + 4*T^2 - 96*T + 120)*x^2 + (8*T^8 - 32*T^7 - 56*T^6 +
      384*T^5 + 8*T^4 - 1952*T^3 + 1240*T^2 + 4800*T - 6000)*x -
      4*T^12 + 16*T^11 + 52*T^10 - 288*T^9 - 272*T^8 + 2528*T^7 +
      104*T^6 - 13056*T^5 + 6620*T^4 + 38800*T^3 - 33500*T^2 - 60000*T
      + 75000, ((-6*T^6 + 48*T^5 - 102*T^4 - 120*T^3 + 768*T^2 -
      1008*T + 360)*x + (12*T^9 - 90*T^8 + 48*T^7 + 1284*T^6 -
      3732*T^5 - 738*T^4 + 17832*T^3 - 24456*T^2 + 1440*T + 12600))*y
      + (3*T^4 - 30*T^3 + 114*T^2 - 180*T + 72)*x^3 + (3*T^8 - 18*T^7
      - 6*T^6 + 228*T^5 - 819*T^4 + 2778*T^3 - 6378*T^2 + 6012*T +
      360)*x^2 + (-12*T^11 + 54*T^10 + 264*T^9 - 1491*T^8 - 1230*T^7 +
      14646*T^6 - 10824*T^5 - 34755*T^4 - 9858*T^3 + 205446*T^2 -
      209340*T - 12600)*x + 12*T^14 - 60*T^13 - 285*T^12 + 1710*T^11 +
      3126*T^10 - 25140*T^9 - 9678*T^8 + 217260*T^7 - 150492*T^6 -
      901140*T^5 + 1169067*T^4 + 2082870*T^3 - 4518750*T^2 + 724500*T
      + 2205000 ];
    u,v,w,x := Explode(ff);
    a := w*x;
    xx := [v - 6*w,(u + T*x)*w,(T - 3)*v - 3*w,T,1];
    x0,x1,x2,x3,x4 := Explode(xx);
    aa := [a, (-2*a*x2*x3 + 2*x1*x2 - 2*a*x0*x4)/(x3^2 - 6*x4^2)];
  when [8,7] : 
    ff := [ ((1152*T^8 + 9216*T^6 - 7776*T^4 - 2592*T^2)*x +
      (9216*T^12 - 44928*T^10 + 65952*T^8 - 28512*T^6 - 3672*T^4 +
      1944*T^2))*y + (-17280*T^10 + 18432*T^8 + 29088*T^6 -
      30240*T^4)*x^2 + (-23040*T^16 + 88704*T^14 - 21888*T^12 -
      288576*T^10 + 440640*T^8 - 240984*T^6 + 45792*T^4 - 648*T^2)*x -
      18432*T^20 + 179712*T^18 - 726912*T^16 + 1572768*T^14 -
      1975104*T^12 + 1468008*T^10 - 635904*T^8 + 155952*T^6 -
      20736*T^4 + 648*T^2, ((-384*T^6 - 3456*T^4 - 864*T^2)*x -
      1536*T^14 + 9216*T^12 - 21888*T^10 + 30336*T^8 - 24192*T^6 +
      8064*T^4)*y + (128*T^12 + 640*T^10 + 2016*T^8 + 3648*T^6 -
      9000*T^4 - 1944*T^2 - 486)*x^2 + (1024*T^18 - 1280*T^16 -
      5120*T^14 + 34752*T^12 - 135744*T^10 + 226512*T^8 - 148800*T^6 +
      35028*T^4 - 7344*T^2 + 972)*x + 2048*T^24 - 21504*T^22 +
      90240*T^20 - 170624*T^18 + 31584*T^16 + 527232*T^14 -
      1165032*T^12 + 1247064*T^10 - 770598*T^8 + 287640*T^6 -
      66204*T^4 + 8640*T^2 - 486, (96*T^10 + 1584*T^8 + 5616*T^6 -
      8424*T^4 - 5346*T^2 - 729)*x^2 + (768*T^16 - 9408*T^14 +
      53856*T^12 - 148896*T^10 + 193536*T^8 - 107676*T^6 + 21222*T^4 -
      4860*T^2 + 1458)*x + 1536*T^22 - 21504*T^20 + 130272*T^18 -
      445392*T^16 + 942480*T^14 - 1279368*T^12 + 1126710*T^10 -
      646137*T^8 + 243648*T^6 - 60750*T^4 + 9234*T^2 - 729 ];
    u,v,x := Explode(ff);
    T0 := 2*T;
    y := (8*T0^3 - 48*T0)/(T0^2 - 6)*x;
    w := -((T0^2 + 12)*u - 3*(T0^2 + 3)*v + (T0^4 - 7*T0^2 + 6)*x)/(18*T0);
    aa := [(1/3)*(-T0^2 + 6)*u*x - 3*u*v - 3*w^2,
           (1/9)*u^2*y + (1/3)*(2*T0^2 - 4)*u*w*x + 2*u*v*w + 2*w^3];
    xx := [(1/3)*T0*u + w,-2*T0*u*v - 2*w^2*T0,2*u + w*T0,T0,1];
  when [9,1] : 
    ff := [ (96*T^7 + 192*T^6 + 168*T^5 + 144*T^4 + 72*T^3)*y - 12*x^3
      + (32*T^5 + 176*T^4 + 144*T^3 + 36*T^2 + 36*T)*x^2 + (192*T^9 +
      96*T^8 - 48*T^7 + 168*T^6 - 24*T^5 - 180*T^4 - 144*T^3 -
      108*T^2)*x + 384*T^13 + 1152*T^12 + 2304*T^11 + 4704*T^10 +
      6696*T^9 + 6336*T^8 + 4752*T^7 + 2592*T^6 + 648*T^5, ((-8*T^3 -
      8*T^2)*x - 24*T^5 - 48*T^4 - 24*T^3)*y + (-2*T + 3)*x^3 +
      (-24*T^5 - 84*T^4 - 62*T^3 - 11*T^2 - 9*T)*x^2 + (144*T^8 +
      96*T^7 - 276*T^6 - 306*T^5 - 93*T^4 + 12*T^3 + 27*T^2)*x -
      96*T^11 - 288*T^10 - 504*T^9 - 960*T^8 - 1296*T^7 - 864*T^6 -
      216*T^5, x^3 + (-12*T^4 - 12*T^3 - 3*T^2 - 3*T)*x^2 + (48*T^8 +
      96*T^7 + 84*T^6 + 72*T^5 + 45*T^4 + 18*T^3 + 9*T^2)*x ];
    T0 := 2*T + 3;
    u,v,w := Explode(ff);
    g0 := u^2 + 3*u*v + 9*u*w + 3*v^2 + 6*T0*v*w + 3*(T0^2 - 12*T0 + 24)*w^2;
    g1 := T0^3 - 9*T0^2 + 36*T0 - 36;
    s := -g0/(4*g1);
    aa := [12*s - 3*w^2, u*s + 2*w^3];
    xx := [6*v*s - 3*w^3,2*s*T0 - w^2,w,1];
  when [9,2] :
    ff := [ ((-2*T - 4)*x + (4*T^7 + 26*T^6 + 78*T^5 + 70*T^4 + 8*T^3
      - 78*T^2 - 74*T - 34))*y + (-14*T^4 - 28*T^3 - 12*T^2 + 2*T -
      2)*x^2 + (4*T^10 + 146*T^9 + 558*T^8 + 720*T^7 + 234*T^6 -
      450*T^5 - 732*T^4 - 348*T^3 - 108*T^2 + 8*T - 32)*x - 24*T^15 +
      456*T^14 + 114*T^13 - 828*T^12 - 2346*T^11 - 816*T^10 + 2436*T^9
      + 4302*T^8 + 1926*T^7 - 2262*T^6 - 3390*T^5 - 1860*T^4 + 648*T^3
      + 978*T^2 + 636*T + 30, ((-T - 2)*x - 6*T^7 - 15*T^6 + 39*T^5 +
      51*T^4 + 60*T^3 - 39*T^2 - 45*T - 45)*y + (-3*T^4 - 6*T^3 -
      6*T^2 - 3*T - 9)*x^2 + (-26*T^10 - 67*T^9 + 108*T^8 + 300*T^7 +
      393*T^6 + 117*T^5 - 162*T^4 - 306*T^3 - 225*T^2 - 112*T - 20)*x
      + 156*T^15 - 12*T^14 - 237*T^13 - 918*T^12 + 30*T^11 + 1200*T^10
      + 1851*T^9 - 18*T^8 - 2097*T^7 - 1686*T^6 - 6*T^5 + 1542*T^4 +
      711*T^3 + 6*T^2 - 408*T - 114, ((-T - 2)*x - 6*T^7 + 9*T^6 +
      39*T^5 + 51*T^4 + 12*T^3 - 39*T^2 - 45*T - 21)*y + (-11*T^4 -
      22*T^3 - 6*T^2 + 5*T + 7)*x^2 + (-58*T^10 + 109*T^9 + 504*T^8 +
      684*T^7 + 153*T^6 - 675*T^5 - 834*T^4 - 354*T^3 + 171*T^2 +
      208*T + 92)*x + 348*T^15 + 372*T^14 - 717*T^13 - 2790*T^12 -
      2478*T^11 + 1608*T^10 + 6015*T^9 + 5202*T^8 - 441*T^7 - 4866*T^6
      - 4458*T^5 - 1074*T^4 + 1107*T^3 + 1362*T^2 + 624*T + 186,
      (8*T^10 + 40*T^9 + 66*T^8 + 24*T^7 - 72*T^6 - 132*T^5 - 72*T^4 +
      24*T^3 + 66*T^2 + 40*T + 8)*x - 48*T^15 - 528*T^14 - 1428*T^13 -
      1728*T^12 + 222*T^11 + 3768*T^10 + 5394*T^9 + 2502*T^8 -
      2736*T^7 - 5334*T^6 - 3558*T^5 - 120*T^4 + 1638*T^3 + 1362*T^2 +
      516*T + 78 ];
    T0 := (2*T - 3)/(2*T + 1);
    u,v,w,s := Explode(ff);
    aa := [3*u*s - 3*w^2, 3*v*s^2 - 3*u*w*s + 2*w^3];
    xx := [-3*s^2*T0 - w^2,s - w,2*w,1];
  when [11,1] : 
    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]>;
    ff := [ (2*T^8 - 39*T^7 + 281*T^6 - 997*T^5 + 2161*T^4 - 2976*T^3
      + 2592*T^2 - 1280*T + 256)*y + (4*T^2 - 16*T + 16)*x^3 + (-2*T^7
      + 21*T^6 - 142*T^5 + 397*T^4 - 554*T^3 + 460*T^2 - 224*T)*x^2 +
      (2*T^11 - 17*T^10 + 114*T^9 - 543*T^8 + 2549*T^7 - 7750*T^6 +
      15087*T^5 - 19430*T^4 + 16432*T^3 - 8416*T^2 + 2304*T - 256)*x -
      108*T^13 + 1200*T^12 - 5935*T^11 + 18109*T^10 - 39412*T^9 +
      65043*T^8 - 83001*T^7 + 81030*T^6 - 58058*T^5 + 28040*T^4 -
      7904*T^3 + 960*T^2, ((4*T^3 - 36*T^2 + 64*T - 64)*x - 8*T^7 +
      108*T^6 - 532*T^5 + 1488*T^4 - 2448*T^3 + 2448*T^2 - 1280*T +
      256)*y + (8*T - 16)*x^3 + (-32*T^5 + 156*T^4 - 544*T^3 + 860*T^2
      - 608*T + 128)*x^2 + (56*T^9 - 436*T^8 + 1868*T^7 - 5424*T^6 +
      11656*T^5 - 16952*T^4 + 15544*T^3 - 8288*T^2 + 2304*T - 256)*x -
      48*T^13 + 592*T^12 - 3372*T^11 + 11832*T^10 - 28624*T^9 +
      51052*T^8 - 69880*T^7 + 73112*T^6 - 55424*T^5 + 27680*T^4 -
      7904*T^3 + 960*T^2, ((-2*T^4 + 22*T^3 - 68*T^2 + 96*T - 64)*x +
      (16*T^7 - 188*T^6 + 718*T^5 - 1610*T^4 + 2280*T^3 - 2096*T^2 +
      1152*T - 256))*y + (-4*T^2 + 16*T - 16)*x^3 + (8*T^6 - 10*T^5 +
      20*T^4 - 138*T^3 + 196*T^2 - 96*T)*x^2 + (-12*T^10 + 54*T^9 -
      104*T^8 - 624*T^7 + 3420*T^6 - 8124*T^5 + 11978*T^4 - 11576*T^3
      + 6944*T^2 - 2176*T + 256)*x + 96*T^13 - 1016*T^12 + 4864*T^11 -
      14614*T^10 + 31660*T^9 - 52384*T^8 + 67294*T^7 - 66604*T^6 +
      49036*T^5 - 24808*T^4 + 7424*T^3 - 960*T^2, (8*T^5 - 104*T^4 +
      448*T^3 - 928*T^2 + 1024*T - 512)*x^2 + (-4*T^10 + 78*T^9 -
      658*T^8 + 3130*T^7 - 8578*T^6 + 14368*T^5 - 14608*T^4 + 8192*T^3
      - 1792*T^2)*x + 24*T^13 - 440*T^12 + 3198*T^11 - 12536*T^10 +
      30078*T^9 - 46758*T^8 + 47212*T^7 - 29578*T^6 + 10176*T^5 -
      1440*T^4 ]; 
    x1,x2,x3,x4 := Explode(ff);
    assert QuadForms(T,x1,x2,x3,x4) eq [0,0];
    a,b,c,d := MobiusCoeffs(x1,x2,x3,x4);
    U := (d*T - b)/(-c*T + a);	 
    x5 := (1/6)*(8*U*x2*x3 - 4*U*x2*x4 + 8*U*x3^2 + 8*U*x3*x4 + 2*x1*x2^2
      + 8*x1*x2*x3 - 3*x2^3 + 13*x2*x3^2 - 3*x2*x3*x4 + 4*x3^3
      + 4*x3^2*x4)/x2^2;
    aa := [ (4*U + 3*x3)*x4 - 3*x5^2,
            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 ];
  when [12,1] :
    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;
    C := Curve(Proj(R),G);
    E1,map1 := EllipticCurve(C,C![1,0,0]);
    flag,map2 := IsIsomorphic(E,E1);
    assert flag;
    mapEtoC := Expand(map2 * map1^(-1));
    xx := DefiningPolynomials(mapEtoC);
    x3 := 2*(x - 4*(T^4 - 1)*(T^2 + 3))^2;
    x1,x2,x3 := Explode([KE!(xx[i]/xx[3])*x3 : i in [1..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 G2*T^2 + G1*T + G0 eq 0;
    T0 := x2/x3;
    x := 2*(T0 + 3)/(T0^2 + 3)*x1/x3;
    Y := (1/2)*(2*G2*T + G1)/x3^3;
    Y *:= (8*T0^3 + 56*T0^2 + 120*T0 + 72)/(T0^2 + 3);
    y := (Y + 4*(T0 - 1)*(T0 + 3)^2*x)/(2*((T0 + 1)^2*x^2 + 4*(T0 + 3)^2));
    y := 8*(T0 + 3)*x*y - (T0 + 3)^2*((T0 - 3)*x^2 + 12*(T0 + 1));
    y /:= 144*(T0 + 1);
    w := (T0 + 3)*(x^2 + 12)/24 - 2*T0*y/(T0 + 3);
    v := 2*(w - 2*y)*y + (T0 + 1)/(T0 - 1)*(x^2*(y + 1) - (w + 2*y)^2);
    a := 2*v + w^2 - x^2;
    aa := [ a, (1/6)*(x^2*y - a)*(x^2 - 12*y) + 8*y^3 ];
    xx := [ x^2 - 6*y, 1, y, v*x, w*x, x];
  end case; 
  return E,aa,xx;  
end function;

// Each point on the surface Z^*(N,eps) in Theorem 3.1 determines
// a pair of elliptic curves. The following function computes those
// two elliptic curves, using the formulae in Section 3.

function CongruentEllipticCurves(N,eps,uu,y)
  assert [N,eps] in {[6,5],[10,1],[10,3]};
  u,v,w := Explode(uu);
  case [N,eps] :
  when [6,5] :
    rr := [(u*w^2 + 4*v^3 - u*(u + 3*v)^2 + e*y)/(8*u) : e in [1,-1]];
    aa := [ -3*r*(r + u*v + 2*v^2) : r in rr];
    bb := [ r*(u + 3*v)*(r*u + v^3) + 2*r^2*(r + 3*v^2) : r in rr ];
    d := u*(u + 3*v);
  when [10,1] : 
    RHS := u*(u^2 - 11*u*v - v^2) + (12*u + v)*w^2;
    rt := 3*u^2 - 4*u*v + 4*w^2;
    ttt := [(RHS + e*y)/(2*w^2): e in [-1,1]];
    aa := [ -3*(8*u - v)*(32*u - v) - 288*rt
          + 30*(28*u + v)*tt - 75*tt^2 : tt in ttt ];
    bb := [ -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 : tt in ttt];
    d := -(35*u^2 - 5*u*v + 4*w^2);
  when [10,3] : 
    r := u^2/(8*u - v + 4*w);
    c := u^2*(11*v + 8*w) + (8*u - v + 4*w)*w^2;
    ss := [(-c + e*y)/(2*(4*u - v + 4*w)^2): e in [-1,1]];
    aa := [ 3*(8*r*s^3 + 4*r*s^2*v + 6*r*s^2*w + r*s*w^2 - s^4) : s in ss];
    bb := [ 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 : s in ss];
    d := r*(2*v - w)*(-4*r*v + 3*r*w + w^2);
  end case;
  E1 := EllipticCurve([aa[1],bb[1]]);
  E2 := EllipticCurve([d^2*aa[2],d^3*bb[2]]);
  return E1,E2;
end function;

// In the cases (N,eps) = (6,5), (10,1), (10,3) we computed Z(N,eps)
// as a degree 3 cover of a K3-surface. This function gives the map
// to the K3-surface.

function MapToZstar(N,eps);
  E := ModularDiagonalQuotientSurface(N,eps);
  KE<x,y> := FunctionField(E);
  case [N,eps] :
  when [6,5] :
    uu := [ -4*x + 8*T^4 + 12*T^3 - 36*T^2 - 32*T + 48, T^2*x - 2*T^6
      - 6*T^5 + 12*T^3 + 12*T^2 - 16, (-T^3 - 3*T^2 + 4)*x + 2*T^7 -
      18*T^5 - 4*T^4 + 36*T^3 - 16*T ];
    y := ((-72*T^7 - 216*T^6 + 96*T^5 + 576*T^4 - 384*T^2)*x +
      (144*T^11 + 576*T^10 - 624*T^9 - 4416*T^8 + 12640*T^6 + 2304*T^5
      - 17280*T^4 - 2048*T^3 + 10752*T^2 - 2048))*y + (-108*T^9 -
      108*T^8 + 792*T^7 + 576*T^6 - 1728*T^5 - 576*T^4 + 1152*T^3)*x^2
      + (144*T^13 - 3216*T^11 - 3456*T^10 + 16704*T^9 + 19296*T^8 -
      40224*T^7 - 38016*T^6 + 51840*T^5 + 29184*T^4 - 32256*T^3 -
      6144*T^2 + 6144*T)*x + 144*T^17 + 1008*T^16 + 1104*T^15 -
      7152*T^14 - 17568*T^13 + 12928*T^12 + 71616*T^11 + 20736*T^10 -
      134720*T^9 - 120384*T^8 + 122112*T^7 + 215040*T^6 - 39936*T^5 -
      202752*T^4 - 4096*T^3 + 98304*T^2 - 16384;
  when [10,1] :
    uu := [ (-2*T^3 + 2*T^2)*x - 80*T^9 + 296*T^8 - 380*T^7 + 168*T^6
      + 20*T^5 - 24*T^4, 8*x - 16*T^7 + 296*T^6 - 780*T^5 + 632*T^4 -
      32*T^3 - 96*T^2, (T^3 + 2*T^2 - 4*T)*x + 32*T^9 - 32*T^8 -
      240*T^7 + 540*T^6 - 340*T^5 - 8*T^4 + 48*T^3 ];
    y := ((-288*T^11 - 48*T^10 + 1960*T^9 - 2640*T^8 + 960*T^7 +
      192*T^6 - 128*T^5)*x + (512*T^18 - 10752*T^17 + 22912*T^16 +
      60544*T^15 - 295520*T^14 + 451840*T^13 - 301536*T^12 +
      39232*T^11 + 55424*T^10 - 22400*T^9 - 1792*T^8 + 1536*T^7))*y +
      (2592*T^14 - 3456*T^13 - 17712*T^12 + 50892*T^11 - 48104*T^10 +
      12592*T^9 + 7104*T^8 - 4032*T^7 - 128*T^6 + 256*T^5)*x^2 +
      (-4608*T^21 + 165888*T^20 - 621952*T^19 - 228352*T^18 +
      5721952*T^17 - 14172192*T^16 + 16359232*T^15 - 8490720*T^14 -
      492128*T^13 + 2723008*T^12 - 889984*T^11 - 198656*T^10 +
      133632*T^9 + 1024*T^8 - 6144*T^7)*x - 110592*T^27 + 2820096*T^26
      - 16211968*T^25 + 26747392*T^24 + 75841536*T^23 - 459305856*T^22
      + 1036935872*T^21 - 1310212352*T^20 + 910279808*T^19 -
      183337088*T^18 - 224659776*T^17 + 181104896*T^16 - 23370240*T^15
      - 27038720*T^14 + 10277376*T^13 + 1215488*T^12 - 1018880*T^11 +
      6144*T^10 + 36864*T^9;
  when [10,3] :
    uu := [ (T^2 + 2*T + 1)*x + 14*T^6 + 46*T^5 + 60*T^4 + 44*T^3 +
      22*T^2 + 6*T, 4*T^3*x + 28*T^7 + 36*T^6 + 48*T^5 + 20*T^4 +
      12*T^3, (T^3 - 2*T^2 - 3*T)*x + 14*T^7 + 18*T^6 - 4*T^5 - 40*T^4
      - 46*T^3 - 26*T^2 - 12*T ];
    y := ((-56*T^10 - 520*T^9 - 1792*T^8 - 2808*T^7 - 1920*T^6 -
      248*T^5 + 480*T^4 + 504*T^3 + 216*T^2)*x - 784*T^14 - 6720*T^13
      - 24704*T^12 - 51440*T^11 - 69088*T^10 - 63776*T^9 - 38336*T^8 -
      9920*T^7 + 6320*T^6 + 9696*T^5 + 6208*T^4 + 2352*T^3 +
      576*T^2)*y + (-28*T^13 - 36*T^12 + 1436*T^11 + 8328*T^10 +
      20416*T^9 + 27360*T^8 + 21104*T^7 + 7128*T^6 - 3076*T^5 -
      5052*T^4 - 2988*T^3 - 864*T^2)*x^2 + (-784*T^17 - 3584*T^16 +
      9232*T^15 + 107856*T^14 + 362144*T^13 + 700848*T^12 +
      911392*T^11 + 833664*T^10 + 521360*T^9 + 172992*T^8 - 48848*T^7
      - 110928*T^6 - 83616*T^5 - 39664*T^4 - 12000*T^3 - 2304*T^2)*x -
      5488*T^21 - 43120*T^20 - 123648*T^19 - 78640*T^18 + 481024*T^17
      + 1776784*T^16 + 3406528*T^15 + 4464336*T^14 + 4273120*T^13 +
      2963184*T^12 + 1315200*T^11 + 95728*T^10 - 428288*T^9 -
      438480*T^8 - 262464*T^7 - 108432*T^6 - 29808*T^5 - 5184*T^4;
  end case;
  return E,uu,y;
end function;
