
// Title : On families of 13-congruent elliptic curves
// Author : T. Fisher
// Version : January 2020

// This is a MAGMA file checking some of the calculations in
// the paper relating to the surfaces Z(13,1) and Z(13,2)

// The calculations are checked in roughly the order of the
// paper, the main exception being that we do not check the
// results stated in Sections 2.2 and 2.3 until later.

MC := MonomialCoefficient;
MD := MonomialsOfDegree;
Deriv := Derivative;

load "congr13-formulae.m";

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

// The group G defined at the start of Section 5.1

K<z13> := CyclotomicField(13);
sqrt13 := &+[LegendreSymbol(a,13)*z13^a : a in [0..12]];
assert sqrt13^2 eq 13;
xi := [z13^i + z13^(-i): i in [1..12]];
I7 := IdentityMatrix(K,7);
M2 := (1/sqrt13)*Matrix(K,7,7,[
    1,   1,     1,     1,     1,     1,     1,
    2, xi[2], xi[4], xi[8], xi[3], xi[6], xi[1],
    2, xi[4], xi[8], xi[3], xi[6], xi[1], xi[2],
    2, xi[8], xi[3], xi[6], xi[1], xi[2], xi[4],
    2, xi[3], xi[6], xi[1], xi[2], xi[4], xi[8],
    2, xi[6], xi[1], xi[2], xi[4], xi[8], xi[3],
    2, xi[1], xi[2], xi[4], xi[8], xi[3], xi[6]
]);
M6 := Matrix(7,7,[-I7[j] : j in [1,7,2,3,4,5,6]]);
M13 := DiagonalMatrix(K,[z13^i : i in [0,1,4,3,12,9,10]]);

G := MatrixGroup<7,K|M2,M13>;
assert #G eq 1092;
assert IsIsomorphic(G,PSL(2,13));
assert M6 in G;

sigma := hom<K->K|z13^2>;
M2tilde := Matrix(7,7,[sigma(M2[i,j]): i,j in [1..7]]);
M13tilde := Matrix(7,7,[sigma(M13[i,j]): i,j in [1..7]]);

// A function for computing a basis for the vector space of
// bi-invariants of degree (d1,d2) (see Definitions 6.1 and 6.4)
// using the efficient averaging method in Remark 5.6.

function BiInvariants(k,RR,d1,d2,dim)
  assert k in [1,2];
  printf "(%o,%o) ",d1,d2;
  R7 := PolynomialRing(K,7);
  mons1 := MD(R7,d1);
  mons2 := MD(R7,d2);
  dd := [Integers(13)|0,1,4,3,12,9,10];
  ww1 := [&+[Degree(m,i)*dd[i] : i in [1..7]]: m in mons1];
  ww2 := [&+[Degree(m,i)*dd[i] : i in [1..7]]: m in mons2];
  P6P6 := ProductProjectiveSpace(K,[6,6]);
  RRK := CoordinateRing(P6P6);
  xx := [RRK.i : i in [1..7]];
  yy := [RRK.i : i in [8..14]];
  mons := [Evaluate(mons1[i],xx)*Evaluate(mons2[j],yy)
    : j in [1..#mons2], i in [1..#mons1] | ww1[i] + k*ww2[j] eq 0];
  M2alt := k eq 1 select M2 else M2tilde;
  M2subst := [Evaluate(R7.i^M2,xx): i in [1..7]]
                  cat [Evaluate(R7.i^M2alt,yy): i in [1..7]];
  M6subst := [Evaluate(R7.i^M6,xx): i in [1..7]]
                  cat [Evaluate(R7.i^M6,yy): i in [1..7]];
  V := VectorSpace(Rationals(),#mons);
  U := sub<V|>;
  ctr := 0;
  vecs := [];
  while Dimension(U) lt dim do
    ff := [Random(mons)];
    for i in [1..5] do
      ff cat:= [Evaluate(ff[#ff],M6subst)];
    end for;
    f1 := &+ff;
    f2 := 13^2*(f1 + 13*Evaluate(f1,M2subst));
    vec := V![MC(f2,m): m in mons];
    printf (vec in U) select "." else "*";
    if vec notin U then
      U +:= sub<V|vec>;
      Append(~vecs,vec);
    end if;
  end while;
  print "";
  return [RR|&+[u[i]*mons[i] : i in [1..#mons]]: u in vecs];
end function;

// Procedure checking a (skew) bi-invariant.

procedure CheckBiInvariant(k,inv:proof:=true);
  assert k in [1,2];
  R7 := PolynomialRing(K,7);
  P6P6 := ProductProjectiveSpace(K,[6,6]);
  RRK := CoordinateRing(P6P6);
  xx := [RRK.i : i in [1..7]];
  yy := [RRK.i : i in [8..14]];
  M2alt := k eq 1 select M2 else M2tilde;
  M13alt := k eq 1 select M13 else M13tilde;
  M2subst := [Evaluate(R7.i^M2,xx) : i in [1..7]];
  M2subst cat:= [Evaluate(R7.i^M2alt,yy) : i in [1..7]];
  M13subst := [Evaluate(R7.i^M13,xx) : i in [1..7]];
  M13subst cat:= [Evaluate(R7.i^M13alt,yy) : i in [1..7]];
  assert Evaluate(RRK!inv,M13subst) eq RRK!inv;
  if proof then
    assert Evaluate(RRK!inv,M2subst) eq RRK!inv;
  end if;
end procedure;

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

// The equations for X(13) in P^6 (see Section 5.1)

R<x0,x1,x2,x3,x4,x5,x6> := PolynomialRing(Rationals(),7);
cubics := [
    -2*x0^3 + x0*(x1*x4 + x2*x5 + x3*x6) + (x1*x3*x5 + x2*x4*x6),
    x0*x1^2 + 2*x0*x3*x4 + 2*x6*x1*x2 + x2*x4^2 + x5*x3^2 + x6*x5^2,
    x0*x2^2 + 2*x0*x4*x5 + 2*x1*x2*x3 + x3*x5^2 + x6*x4^2 + x1*x6^2,
    x0*x3^2 + 2*x0*x5*x6 + 2*x2*x3*x4 + x4*x6^2 + x1*x5^2 + x2*x1^2,
    x0*x4^2 + 2*x0*x6*x1 + 2*x3*x4*x5 + x5*x1^2 + x2*x6^2 + x3*x2^2,
    x0*x5^2 + 2*x0*x1*x2 + 2*x4*x5*x6 + x6*x2^2 + x3*x1^2 + x4*x3^2,
    x0*x6^2 + 2*x0*x2*x3 + 2*x5*x6*x1 + x1*x3^2 + x4*x2^2 + x5*x4^2
];
quartics := [ x0^4 - x0^2*x1*x4 + x1*x3*x4*x6, x0^4 - x0^2*x2*x5 +
    x1*x2*x4*x5, x0^4 - x0^2*x3*x6 + x2*x3*x5*x6, x0^3*x1 + x0*x3*x4^2
    + x3^2*x4*x5, x0^3*x2 + x0*x4*x5^2 + x4^2*x5*x6, x0^3*x3 +
    x0*x5*x6^2 + x1*x5^2*x6, x0^3*x4 + x0*x1^2*x6 + x1*x2*x6^2,
    x0^3*x5 + x0*x1*x2^2 + x1^2*x2*x3, x0^3*x6 + x0*x2*x3^2 +
    x2^2*x3*x4, x0*x1^2*x5 + x0*x3*x4*x5 - x1*x2*x3^2, x0*x2^2*x6 +
    x0*x4*x5*x6 - x2*x3*x4^2, x0*x3^2*x1 + x0*x1*x5*x6 - x3*x4*x5^2,
    x0*x4^2*x2 + x0*x1*x2*x6 - x4*x5*x6^2, x0*x5^2*x3 + x0*x1*x2*x3 -
    x1^2*x5*x6, x0*x6^2*x4 + x0*x2*x3*x4 - x1*x2^2*x6, x0^3*x1 +
    x0^2*x2*x6 - x2^2*x5*x6 - x4*x5^2*x6, x0^3*x2 + x0^2*x1*x3 -
    x1*x3^2*x6 - x1*x5*x6^2, x0^3*x3 + x0^2*x2*x4 - x1^2*x2*x6 -
    x1*x2*x4^2, x0^3*x4 + x0^2*x3*x5 - x1*x2^2*x3 - x2*x3*x5^2,
    x0^3*x5 + x0^2*x4*x6 - x2*x3^2*x4 - x3*x4*x6^2, x0^3*x6 +
    x0^2*x1*x5 - x1^2*x4*x5 - x3*x4^2*x5, x0^2*x1^2 + x0^2*x3*x4 +
    x0*x3^2*x5 - x1^2*x3*x6, x0^2*x2^2 + x0^2*x4*x5 + x0*x4^2*x6 -
    x1*x2^2*x4, x0^2*x3^2 + x0^2*x5*x6 + x0*x1*x5^2 - x2*x3^2*x5,
    x0^2*x1*x6 + x0^2*x4^2 + x0*x2*x6^2 - x3*x4^2*x6, x0^2*x1*x2 +
    x0^2*x5^2 + x0*x1^2*x3 - x1*x4*x5^2, x0^2*x2*x3 + x0^2*x6^2 +
    x0*x2^2*x4 - x2*x5*x6^2, x0^2*x1*x2 + x0^2*x5^2 + x0*x4*x5*x6 -
    x1*x2*x3*x6, x0^2*x2*x3 + x0^2*x6^2 + x0*x1*x5*x6 - x1*x2*x3*x4,
    x0^2*x1^2 + x0^2*x3*x4 + x0*x1*x2*x6 - x2*x3*x4*x5, x0^2*x2^2 +
    x0^2*x4*x5 + x0*x1*x2*x3 - x3*x4*x5*x6, x0^2*x3^2 + x0^2*x5*x6 +
    x0*x2*x3*x4 - x1*x4*x5*x6, x0^2*x1*x6 + x0^2*x4^2 + x0*x3*x4*x5 -
    x1*x2*x5*x6, x0^2*x1*x4 + x1*x2^3 + x2^2*x5^2 + x4*x5^3,
    x0^2*x2*x5 + x2*x3^3 + x3^2*x6^2 + x5*x6^3, x0^2*x3*x6 + x1^3*x6 +
    x1^2*x4^2 + x3*x4^3, x0^2*x1*x3 + x1*x2*x3*x5 + x1*x5*x6^2 +
    x3*x5^3, x0^2*x2*x4 + x1^2*x2*x6 + x2*x3*x4*x6 + x4*x6^3,
    x0^2*x3*x5 + x1^3*x5 + x1*x2^2*x3 + x1*x3*x4*x5, x0^2*x4*x6 +
    x2^3*x6 + x2*x3^2*x4 + x2*x4*x5*x6, x0^2*x1*x5 + x1*x3^3 +
    x1*x3*x5*x6 + x3*x4^2*x5, x0^2*x2*x6 + x1*x2*x4*x6 + x2*x4^3 +
    x4*x5^2*x6 ];

assert #cubics eq 7;
assert #quartics eq 42;
CC := Curve(Proj(R),cubics cat quartics);
assert Degree(CC) eq 42;

// Some invariants and covariants defined in Section 5

Q := x0^2 + x1*x4 + x2*x5 + x3*x6;
F := 2*x0^4 + 6*x0*(x1*x3*x5 + x2*x4*x6)
  + 3*(x1*x2*x4*x5 + x1*x3*x4*x6 + x2*x3*x5*x6)
  + x1*x2^3 + x2*x3^3 + x3*x4^3 + x4*x5^3 + x5*x6^3 + x6*x1^3;
N := Matrix(R,7,7,[
    [ 0, x4, -x5, x6, -x1, x2, -x3 ],
    [ -x4, 0, 0, x5, x0, -x3, 0 ],
    [ x5, 0, 0, 0, -x6, -x0, x4 ],
    [ -x6, -x5, 0, 0, 0, x1, x0 ],
    [ x1, -x0, x6, 0, 0, 0, -x2 ],
    [ -x2, x3, x0, -x1, 0, 0, 0 ],
    [ x3, 0, -x4, -x0, x2, 0, 0 ]
]);

H := func<f|Matrix(7,7,[Deriv(Deriv(f,i),j): i,j in [1..7]])>;
HQinv := H(Q)^(-1);
Compose := func<v,w|Vector(7,[Evaluate(f,Eltseq(w)): f in Eltseq(v)])>;
Dot := func<v,w|&+[H(Q)[i,j]*v[i]*w[j] : i,j in [1..7]]>;
nablaQ := func<f|Vector([Deriv(f,i) : i in [1..7]])*HQinv>;

v3 := nablaQ(F);
v4 := v3*(-N)*HQinv;
v9 := Compose(v3,v3);
c6 := Dot(v4,v9);

function GetCovariant(f0,f1)
  v := [f0,f1];
  M := ChangeRing(M6^(-1),Rationals());
  for i in [1..5] do
    Append(~v,-v[#v]^M);
  end for;
  return Vector(v);
end function;

g0 := 4*x0*(x1*x3*x5 - x2*x4*x6) + x1*x2^3 - x2*x3^3 + x3*x4^3
   - x4*x5^3 + x5*x6^3 - x6*x1^3;
g1 := 4*x0^2*x1^2 - 4*x0^2*x3*x4 + 4*x0*x1*x2*x6 - 2*x0*x3^2*x5
   - 2*x0*x5^2*x6 - x1^3*x4 + x1^2*x2*x5 - 2*x1*x3*x4^2 - x1*x5^3
   - x2^3*x3 - 2*x2^2*x6^2 + 4*x2*x3*x4*x5 + x3^2*x4*x6 + x4*x5*x6^2;
w4 := GetCovariant(g0,g1);

printf "\nComputing the points above j = 0, 1728, infty\n\n";

// We compute the G-orbits of the points recorded in the proof
// of Theorem 5.8, and in the paragraph following the proof.

P<x> := PolynomialRing(Rationals());
n := func<v|v/v[i] where i is Maximum([j : j in [1..7] | v[j] ne 0])>;
rep_lists := <>;

for ctr in [1..3] do
  case ctr :
  when 1 :
    print "Above j = infty";
    pt13 := [ 0, 1, 0, 0, 0, 0, 0 ];
    L<z13> := K;
  when 2 :
    print "Above j = 0";
    L4<a> := NumberField(x^4 + 5*x^2 + 13);
    w := (a^2 + 1)/3;
    assert w^2 + w + 1 eq 0;
    assert a^2 eq -1 + 3*w;
    pt3 := [ -2, a + w, -a + w, a + w, -a + w, a + w, -a + w ];
    L<z13> := ext<L4|CyclotomicPolynomial(13)>;
    Embed(K,L,z13);
  when 3 :
    print "Above j = 1728";
    L6<b> := NumberField(x^6 - 2*x^5 + 2*x^3 - x^2 + 1);
    i := b^4 - 2*b^3 + b^2 + b - 1;
    assert i^2 + 1 eq 0;
    assert b^3 - (i + 1)*b^2 - b + i eq 0;
    pt2 := [ 1, b, b^2 - b, -b^2 + i + 1, b, b^2 - b, -b^2 + i + 1 ];
    L<z13> := ext<L6|CyclotomicPolynomial(13)>;
    Embed(K,L,z13);
  end case;
  sigma := hom<L->L|z13^2>;
  MM2 := ChangeRing(Transpose(M2),L);
  MM6 := ChangeRing(Transpose(M6),L);
  MM13 := ChangeRing(Transpose(M13),L);
  case ctr :
  when 1 : 
    pts := {n(Vector(L,pt13))};
    pts := {n(P*MM2^r): r in [0..1],P in pts};
    pts := {n(P*MM6^r): r in [0..5],P in pts};
    pts := {n(P*MM13^r): r in [0..12],P in pts};
    assert #pts eq 84;
  when 2 : 
    pts := {n(Vector(L,pt3))};
    pts := {n(P*MM13^r): r in [0..12],P in pts};
    pts := {n(P*MM2^r): r in [0..1],P in pts};
    pts := {n(P*MM6^r): r in [0..1],P in pts};
    pts := {n(P*MM13^r): r in [0..12],P in pts};
    assert #pts eq 364;
  when 3 : 
    pts := {n(Vector(L,pt2))};
    pts := {n(P*MM13^r): r in [0..12],P in pts};
    pts := {n(P*MM2^r): r in [0..1],P in pts};
    pts := {n(P*MM6^r): r in [0..2],P in pts};
    pts := {n(P*MM13^r): r in [0..12],P in pts};
    assert #pts eq 546;
  end case;
  printf "   We have %o points\n",#pts;
  fld := ["Q","L4","L6"][ctr];
  printf "   Dividing into Gal(%o(z13)/%o)-orbits\n",fld,fld;
  printf "   Orbit sizes : ";
  reps := [];
  pts_seen := {};
  for P in pts do
    if P notin pts_seen then
      PP := [];
      P1 := Eltseq(P);
      for loopvar in [1..12] do
        P1 := [sigma(x): x in P1];
        Append(~PP,n(Vector(P1)));
      end for;
      _,m := Minimum([#Sprint(v): v in PP]);
      Append(~reps,Eltseq(PP[m]));
      pts_seen join:= Set(PP);
      printf "%o ",#Set(PP);
    end if;
  end for;
  printf "\n   #orbits = %o\n\n",#reps;
  Append(~rep_lists,reps);
end for;

// We compute a list of pairs of points above j = 0,1728,infty.
// Our list includes enough pairs to be sure that any
// symmetric (or anti-symmetric) bi-invariant defined over Q
// vanishing at these pairs must in fact vanish at all pairs.

pairs := <P1 cat P2 : P2 in rep_lists[1], P1 in <pt13,pt3,pt2>>;
pairs cat:= <pt3 cat P2 : P2 in rep_lists[2]>;
pairs cat:= <pt2 cat P2 : P2 in rep_lists[3]>; 
PP<T> := PolynomialRing(L4);
L24<b> := ext<L4|T^6 - 2*T^5 + 2*T^3 - T^2 + 1>;
Embed(L6,L24,b);
L<z13> := ext<L24|CyclotomicPolynomial(13)>;
Embed(Universe(rep_lists[2,1]),L,z13);
pairs cat:= <pt2 cat P2 : P2 in rep_lists[2]>;
assert #pairs eq 154;

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

print "Section 6\n";

P6P6 := ProductProjectiveSpace(Rationals(),[6,6]);
RR<x0,x1,x2,x3,x4,x5,x6,y0,y1,y2,y3,y4,y5,y6> := CoordinateRing(P6P6);
xx := [x0,x1,x2,x3,x4,x5,x6];
yy := [y0,y1,y2,y3,y4,y5,y6];
xy := Eltseq(Vector(xx) + Vector(yy));
iota1 := yy cat xx;
iota2 := yy cat [-x0,-x2,-x3,-x4,-x5,-x6,-x1];

c6x := Evaluate(c6,xx);
c6y := Evaluate(c6,yy);
eqns := [Evaluate(f,xx): f in cubics cat quartics];
eqns cat:= [Evaluate(f,yy): f in cubics cat quartics];
II := ideal<RR|eqns>;

function RandomPoint(C)
  k := BaseField(C);
  R := CoordinateRing(Ambient(C));
  repeat
    repeat 
      linearform := &+[Random(k)*R.i: i in [1..Rank(R)]];
    until linearform ne 0;
    X := C meet Scheme(Ambient(C),linearform);
    pts := Points(X);
  until #pts gt 0;
  return C!Random(pts);
end function;

p := NextPrime(10^3);
k := GF(p);

R<[A]> := PolynomialRing(k,7);
eqns := [R!f : f in cubics cat quartics];
C := Curve(Proj(R),eqns);
assert Degree(C) eq 42;

printf "Computing some random points over GF(%o)\n",p;
ptsC := [RandomPoint(C): i in [1..30]];

print "";
print "Section 6.1 : THE DIRECT CASE";
print "";

// The Q_ij and F_ij defined at the start of Section 6.1

Q20 := Evaluate(Q,xx);
Q11 := &+[Deriv(Q20,RR.i)*RR.(i+7): i in [1..7]];
Q02 := (1/2)*&+[Deriv(Q11,RR.i)*RR.(i+7): i in [1..7]];
assert Q02 eq Evaluate(Q20,iota1);
assert Q11 eq Evaluate(Q11,iota1);
assert Evaluate(Q,xy) eq Q20 + Q11 + Q02;

F40 := Evaluate(F,xx);
F31 := &+[Deriv(F40,RR.i)*RR.(i+7): i in [1..7]];
F22 := (1/2)*&+[Deriv(F31,RR.i)*RR.(i+7): i in [1..7]];
F13 := (1/3)*&+[Deriv(F22,RR.i)*RR.(i+7): i in [1..7]];
F04 := (1/4)*&+[Deriv(F13,RR.i)*RR.(i+7): i in [1..7]];
assert F04 eq Evaluate(F40,iota1);
assert F13 eq Evaluate(F31,iota1);
assert F22 eq Evaluate(F22,iota1);
assert Evaluate(F,xy) eq F40 + F31 + F22 + F13 + F04;

for inv in [Q20,Q11,Q02,F31,F13] do
  CheckBiInvariant(1,inv);
end for;

// A basis for the space of bi-invariants of degree (3,3)

invs33 := BiInvariants(1,RR,3,3,10);
mons := [ x2^3*y0^2*y1, x1*x4^2*y0^2*y1, x2^2*x3*y0*y1^2,
    x3*x4*x5*y0*y1^2, x2*x3^2*y1^3, x4^3*y1^3, x0^2*x6*y1^3,
    x1*x4*x6*y1^3, x2*x5*x6*y1^3, x3*x6^2*y1^3 ];
coeffmat := Matrix(Rationals(),10,10,[ 1, -2, -2, 4, 0, 1, 1, 1, -1,
  0, -2, 0, 4, -2, -1, -1, -2, -1, 1, 0, 0, -3, 0, 1, -1, 0, 0, 0, 0,
  0, 0, 0, -2, 0, 0, 0, 0, 0, -1, 0, 1, 0, -3, 0, 0, 0, 1, 0, -1, 0,
  1, -1, -3, 0, 0, 0, 1, 0, -1, 0, 1, -1, -2, 4, 0, 0, 1, 1, -1, 0, 0,
  -4, 2, 0, -1, 0, 0, 0, 1, -1, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, -4,
  0, 4, -1, 0, 1, 1, 0, 0 ]);
mat := Matrix([[MC(f,mon): mon in mons]: f in invs33]);
BB := coeffmat*mat^(-1);
invs33 := [&+[BB[i,j]*invs33[j] : j in [1..10]]: i in [1..10]];

// Quadratic relations and the rational surface Sigma

R<z1,z2,z3,z4,z5,z6,z7,z8,z9> := PolynomialRing(Rationals(),10);
LHS := [Q11^3, Q11*Q20*Q02, Q11*F22, Q20*F13 + Q02*F31];
RHS := [
    z1 + z5 - z6 - z7,
    z5 - z6,
    -3*(z6 - z7 + z9),
    z1 + z2 - 4*z4 + z6 + z7 - z8 - 3*z9
];
assert forall{i : i in [1..4] | LHS[i] eq Evaluate(RHS[i],invs33)};

quads := [
  z1*z4 - z3*z5,
  z1*z6 - z3*z7,
  z4*z7 - z5*z6,
  z1^2 + z1*z2 - z2*z4,
  z1*z7 - z1*z8 + z2*z7 - z4*z6 + z5*z7,
  z1*(z1 + z2 - z3 - z4 + z5 - z7) - z2*z6 + z3*z6,
  z4*(z1 + z2 - z3 - z4 + z5) - z3*z8,
  z5*(z1 + z2 - z3 - z4 + z5) - z1*z8,
  z8*(z1 + z5 - z6 - z7) - z5*z7
];

procedure CheckRelations(relns,invs:proof := "heuristic")
  print "Proof method : ",proof;
  case proof :
  when "groebner" :
    for reln in relns do
      time assert Evaluate(reln,invs) in II;
    end for;
  when "heuristic" :
    for ctr in [1..30] do
      P1 := Eltseq(Random(ptsC));
      P2 := Eltseq(Random(ptsC));
      ii := [Evaluate(f,P1 cat P2): f in invs];
      assert forall{r : r in relns | Evaluate(r,ii) eq 0};
      printf ".";
    end for;
    print "";
  when "bezout" :
    for subst in pairs do
      d := AbsoluteDegree(Universe(subst));
      ii := [Evaluate(f,subst) : f in invs];
      assert forall{r : r in relns | Evaluate(r,ii) eq 0};
      printf "%o",[".",".",".","*"][Position([12,48,72,288],d)];
    end for;
    print "";    
  end case;
end procedure;

print "Checking relations to get the surface";

// time CheckRelations([z9] cat quads,invs33:proof := "groebner");
time CheckRelations([z9] cat quads,invs33:proof := "heuristic");
// time CheckRelations([z9] cat quads,invs33:proof := "bezout");

K<r,s> := FunctionField(Rationals(),2);
param := [ r, s, 1, r*(r + s)/s, r^2*(r + s)/s,
     r*f/(s*(r^2 + s - 1)), r^2*f/(s*(r^2 + s - 1)),
     r*(r + s)*f/s^2, 0, 0 ]
  where f is r^3 + r^2*s - r^2 + s^2 - s;
assert forall{q : q in quads | Evaluate(q,param) eq 0};

// A further relation used to compute Z(13,1) as a double
// cover of the rational surface Sigma

SS<u,w,z1,z2,z3,z4,z5,z6,z7,z8> := PolynomialRing(Rationals(),10);
g := z1*z2 + z2^2 - 48*z2*z3 + 48*z2*z5 + 126*z2*z6 - 48*z2*z7 -
    2*z2*z8 + 48*z3^2 - 7*z3*z4 - 57*z3*z5 - 108*z3*z6 + 126*z3*z7 +
    21*z3*z8 + z4^2 - 41*z4*z5 - 26*z4*z6 + 8*z4*z8 + 104*z5^2 -
    60*z5*z6 - 106*z5*z7 + 120*z5*z8 + 37*z6^2 - 10*z6*z7 - 158*z6*z8
    + z7^2 - 70*z7*z8 + z8^2;
invs := [Q20*Q02,Q20*F13 - Q02*F31] cat [invs33[i] : i in [1..8]];
time CheckRelations([w^2 + 64*u^3 - g],invs:proof := "heuristic");
// time CheckRelations([w^2 + 64*u^3 - g],invs:proof := "bezout");

g1,h1 := Z13Equations(1,r,s);
F1 := h1^2 + 4*g1;
pp := Evaluate(RHS[1],param); // Q11^3
qq := Evaluate(RHS[2],param); // Q11*Q20*Q02
gg := Evaluate(g,[0,0] cat [param[i] : i in [1..8]]);
assert gg - 64*qq^3/pp eq ((r - 1)/(s^2*(r^2 + s - 1)))^2*F1;

// Some more bi-invariants and some more relations

print "Checking relations to get the j-maps";

vv := v4*H(Q);
I41 := &+[Evaluate(vv[i],xx)*RR.(i+7): i in [1..7]];
I32 := (-1/3)*&+[Deriv(I41,RR.i)*RR.(i+7): i in [1..7]];
I23 := Evaluate(I32,iota1);

invs := [Q20,Q11,Q02,F31,F13,I32,I23]
         cat [invs33[i] : i in [1..8]] cat [c6x,c6y];

SS<Q20,Q11,Q02,F31,F13,I32,I23,z1,z2,z3,z4,z5,z6,z7,z8,c6x,c6y>
            := PolynomialRing(Rationals(),17);
z := [z1,z2,z3,z4,z5,z6,z7,z8];
alpha := Q20*I23^2;
alpha1 := Q02*I32^2;
beta := Q02^6*I23*c6x; 
beta1 := Q20^6*I32*c6y;
gamma := (Q20*Q02)^3;
f1,f2,f3,f7,g2,h2,l2,l4,l6 := ExtraRelationsDirect(z);
relns := [  f1*(alpha + alpha1) - f3,
            f2*(beta + beta1) - f7,
            Q11*I23*I32 - g2,
            h2*c6x*c6y - Q11*(l6 + l4*gamma + l2*gamma^2 - 64*gamma^3) ];
time CheckRelations(relns,invs:proof := "heuristic");
// time CheckRelations(relns,invs:proof := "bezout");

// Deducing the formulae for j and j' (see Theorem 2.4(iii))
// from the relations above and the parametrisation of Sigma

P<X> := PolynomialRing(K);
L<c> := quo<P|X^2 - F1>;
param := [0 : i in [1..7]] cat param;

asum := Evaluate(f3/f1,param);
aprod := (qq/pp)*Evaluate(g2,param)^2;
adisc := asum^2 - 4*aprod;
flag,rho := IsSquare(adisc/F1);
assert flag;
alpha := (1/2)*(asum + c*rho);
assert Trace(alpha) eq asum;
assert Norm(alpha) eq aprod;

bsum := Evaluate(f7/f2,param);
bprod := Evaluate(l6,param) + Evaluate(l4,param)*gamma
  + Evaluate(l2,param)*gamma^2 - 64*gamma^3 where gamma is qq^3/pp;
bprod *:= Evaluate(g2/h2,param)*(qq^3/pp)^2;
bdisc := bsum^2 - 4*bprod;
flag,rho := IsSquare(bdisc/F1);
assert flag;
beta := (1/2)*(bsum + c*rho); // beware the sign choice
assert Trace(beta) eq bsum;
assert Norm(beta) eq bprod;

j := 1728 - beta^2/(alpha*gamma^4) where gamma is qq^3/pp;
// checking the answer agrees with that in congr13-formulae.m
A1,B1,D1 := Z13ModuliInterpretation(1,r,s);
assert Norm(j) eq A1^3/D1;
assert Norm(j - 1728) eq B1^2/D1;

print "";
print "Section 6.2 : THE SKEW CASE";
print "";

Q20 := Evaluate(Q,xx);
Q02 := Evaluate(Q,yy);

// A basis for the space of skew bi-invariants of degree (2,2)

invs22 := BiInvariants(2,RR,2,2,3);
mons := [x0^2*y0^2,x1*x4*y0^2,x1*x6*y0*y1];
coeffmat := Matrix(Rationals(),3,3,[ 2, 0, 0, 1, -1, 0, -5, 1, -2 ]);
mat := Matrix([[MC(f,mon): mon in mons]: f in invs22]);
BB := coeffmat*mat^(-1);
invs22 := [&+[BB[i,j]*invs22[j] : j in [1..3]]: i in [1..3]];
assert forall{f : f in invs22 | Evaluate(f,iota2) eq f};

for inv in [Q20,Q02] cat invs22 do
  CheckBiInvariant(2,inv);
end for;

// A basis for the space of skew bi-invariants of degree (3,3)

invs33 := BiInvariants(2,RR,3,3,9);
mons := [ x2^2*x3*y0^2*y1, x1*x2^2*y0*y1^2, x3^3*y0*y1^2,
    x0^2*x5*y0*y1^2, x1*x3^2*y1^3, x2^2*x4*y1^3, x4^2*x5*y1^3,
    x1*x5*x6*y1^3, x0*x6^2*y1^3 ];
coeffmat := Matrix(Integers(),9,9,[ 0, -2, 0, 0, -1, 0, 0, -1, 0,
  -3, 1, 0, -1, 0, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 1, 2, 0, -1, 1, 0, 1,
  0, 0, 0, 0, 0, -1, -2, -1, -3, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -3, -1, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, -1, -1, -1, -2, -1 ]);
mat := Matrix([[MC(f,mon): mon in mons]: f in invs33]);
BB := coeffmat*mat^(-1);
invs33 := [&+[BB[i,j]*invs33[j] : j in [1..9]]: i in [1..9]];

// Some relations of small degree, and the rational surface Sigma

print "Checking relations to get the surface";

SS<[u]> := PolynomialRing(Rationals(),9);
time CheckRelations([u[6],u[7],u[8],u[9]],invs33:proof := "heuristic");
// time CheckRelations([u[6],u[7],u[8],u[9]],invs33:proof := "bezout");

SS<t1,t2,t3,u1,u2,u3,u4,u5> := PolynomialRing(Rationals(),[2,2,2,3,3,3,3,3]);
relns := [ u3*u4 - u1*u5,
  u2^3 - u1*u2*u3 - 2*u2^2*u3 - u2^2*u4 + u1*u3*u4 + u2*u3*u4 - u2*u3*u5,
  t2*u4 - t3*u1,
  t1*u1*(u2 - u3) - t2*(u1*u2 - u1*u4 + u2*u3 + u2*u5),
  t1*t2^2 - u1*u3 ];
invs := invs22 cat [invs33[i] : i in [1..5]];
time CheckRelations(relns,invs:proof := "heuristic");
// time CheckRelations(relns,invs:proof := "bezout");

K1<tau> := quo<P|X^3 - (r^2 + r*s + r + 1)*(r^2*s + r*s^2 + r*s + 1)^2>;
param := [-s*(r^2 + r*s + r + 1)/tau,
           r*(r^2*s + r*s^2 + r*s + 1)/tau,
	   (r^2*s + 2*r*s + 1)/tau,
	   r,1,-r*s,f,-s*f]
	   where f is (r^2*s + 2*r*s + 1)/(r^2*s + r*s^2 + r*s + 1);
assert forall{ i : i in [1..5] | Evaluate(relns[i],param) eq 0 };

// Some more bi-invariants and some more relations

print "Checking relations to get the j-maps";

// I32 = 12*(x1*x3*x5 - x2*x4*x6)*y0^2 + .....
I32 := BiInvariants(2,RR,3,2,1)[1];
I32 *:= 12/MC(I32,x1*x3*x5*y0^2);
I23 := Evaluate(I32,iota2);

ff := [Evaluate(f,xx): f in Eltseq(w4*H(Q))];
I41 := &+[RR.(i+7)*ff[i] : i in [1..7]];
I14 := Evaluate(I41,iota2);

c6x := Evaluate(c6,xx);
c6y := Evaluate(c6,yy);

invs := [Q20,Q02,I41,I14] cat invs22
     cat [invs33[i] : i in [1..5]] cat [I32*I23,c6x,c6y];
SS<Q20,Q02,I41,I14,t1,t2,t3,u1,u2,u3,u4,u5,I55,c6x,c6y>
     := PolynomialRing(Rationals(),15);
alpha := Q20^3*I14^2;
alpha1 := Q02^3*I41^2;
beta := Q02^5*I14*c6x;
beta1 := Q20^5*I41*c6y;
gamma := I55;
f5,f8,f10,g10,f11,f12,f13,f15,f18,f20,f23
  := ExtraRelationsSkew([t1,t2,t3,u1,u2,u3,u4,u5]);
relns := [
  Q20*Q02 - (t1 - t2),
  I41*I14 - (t1 - t2)*(u1 + u2 - 3*u3 - u4),
  u1*(u1 + u3)*gamma - f11,
  (t1 + t3)*t3*(alpha + alpha1) - f12,
  4*(t1^2*t3 + u4*u5)*(beta + beta1)
    - (f20 + f15*gamma + f10*gamma^2 + f5*gamma^3),
  g10*c6x*c6y - (f23 + f18*gamma + f13*gamma^2 + f8*gamma^3) ];
time CheckRelations(relns,invs:proof := "heuristic");
// time CheckRelations(relns,invs:proof := "bezout");

// Deducing the formulae for j and j' (see Theorem 2.4(iii))
// from the relations above and the parametrisation of Sigma

g2,h2 := Z13Equations(2,r,s);
F2 := h2^2 + 4*g2;
L<c> := quo<P|X^2 - F2>;
param := [0,0,0,0] cat param cat [0,0,0];
param[13] := Evaluate(f11/(u1*(u1 + u3)),param);

asum := K!(tau*Evaluate(f12/((t1 + t3)*t3),param));
xi := Evaluate((t1 - t2)*(u1 + u2 - 3*u3 - u4),param);
mu := Evaluate(t1 - t2,param);
aprod := K!(tau^2*mu^3*xi^2);
adisc := asum^2 - 4*aprod;
flag,rho := IsSquare(adisc/F2);
assert flag;
alpha_tau := (1/2)*(asum + c*rho);
assert Trace(alpha_tau) eq asum;
assert Norm(alpha_tau) eq aprod;

num := f20 + f15*I55 + f10*I55^2 + f5*I55^3;
den := 4*(t1^2*t3 + u4*u5);
bsum := K!(tau*Evaluate(num/den,param));
num := f23 + f18*I55 + f13*I55^2 + f8*I55^3;
eta := Evaluate(num/g10,param);
bprod := K!(tau^2*mu^5*xi*eta);
bdisc := bsum^2 - 4*bprod;
flag,rho := IsSquare(bdisc/F2);
assert flag;
beta_tau := (1/2)*(bsum - c*rho); // beware the sign choice
assert Trace(beta_tau) eq bsum;
assert Norm(beta_tau) eq bprod;

j := 1728 - beta_tau^2/(alpha_tau*K!(tau*mu^10));
// checking the answer agrees with that in congr13-formulae.m
A2,B2,D2 := Z13ModuliInterpretation(2,r,s);
assert Norm(j) eq A2^3/D2;
assert Norm(j - 1728) eq B2^2/D2;

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

print "\nSection 7\n";

// The formulae for Z(13,k) and the j-maps from Theorem 2.4

R<x,y,z> := PolynomialRing(Rationals(),3);
P2 := Proj(R);
g1,h1 := Z13Equations(1,x/z,y/z);
F1 := R!(z^12*(4*g1 + h1^2));
g2,h2 := Z13Equations(2,x/z,y/z);
F2 := R!(z^14*(4*g2 + h2^2));
F := [F1,F2];

A1,B1,D1 := Z13ModuliInterpretation(1,x/z,y/z);
jprod1 := A1^3/D1;
jsum1 := (jprod1 + 1728^2 - B1^2/D1)/1728;
A2,B2,D2 := Z13ModuliInterpretation(2,x/z,y/z);
jprod2 := A2^3/D2;
jsum2 := (jprod2 + 1728^2 - B2^2/D2)/1728;
jsum := [jsum1,jsum2];
jprod := [jprod1,jprod2];

Kt<t> := FunctionField(Rationals());
P<e> := PowerSeriesRing(Kt);

// Table 1 : Copies of X_0(m) on Z(13,1) and Z(13,2)

Table1 := [
<2, 0, [4, -3 - e^2 + t*e^4, -2 + 2*e], 2^18*(4*t + 1), 4>,
<3, 0, [1, 1 + e + t*e^3, 2 + e], 16*(54*t + 1), 4>,
<4, 0, [1, t*e^2, -1 + e], 4*(32*t + 1), 2>,
<5, 0, [t*e^2, 1, -1 + e], -(16*t^2 + 44*t - 1), 4>,
<6, 0, [e^3, t ,e^2], t^6*(t^2 + 34*t + 1), 18>,
<7, 0, [1, -1 + t*e^2 + t*e^3, -1 + t*e^2], t^4*(t + 1)*(t - 27), 12>,
<8, 0, [e^2 - e^3 - t*e^4, -1, -e], t^2 + 28*t + 68, 14>,
<9, 0, [e, 1 - e^2 + t*e^3, 1], t^2 - 18*t - 27, 6>,
<10, 0, [0, t, 1], t^4*(t - 1)^2*(t^2 + 16*t - 16), 0>,
<11, 1, [1, -e^3 + e^4 + t*e^5, e], (t + 2)*(t^3 - 14*t^2 - 12*t - 4), 20>,
<12, 0, [1, -1 + t*e, e], t^2*(t^2 + 14*t + 1), 4>,
<14, 1, [e, 1, -e^2 + t*e^3], t^4 + 14*t^3 + 19*t^2 + 14*t + 1, 12>,
<15, 1, [1, -1 + t*e, e], (t^2 + t - 1)*(t^2 + 13*t + 11), 4>,
<16, 0, [1, e, 1 + t*e], (4*t^2 + 4*t - 7), 2>,
<17, 1, [t*e, e^2, t], t^8*(t^4 - 6*t^3 - 27*t^2 - 28*t - 16), 8 >,
<18, 0, [t, 0, 1], t^2 + 10*t + 1, 0>,
<19, 1, [-1, t, 1], (t - 1)^4*(t - 2)*(t^3 + 10*t^2 + 12*t + 4), 0>,
<20, 1, [t, 1, 1], (t + 1)^6*(t^4 + 12*t^3 + 28*t^2 + 32*t + 16), 0>,
<21, 1, [1, -t^2, t], t^8*(t + 1)^4*(t^4 + 6*t^3 - 17*t^2 + 6*t + 1), 0>,
<24, 1, [t^3 + 2*t^2, -1, t^2 + 2*t],
  t^10*(t + 1)^6*(t + 2)^6*(t^4 - 4*t^3 - 16*t^2 - 8*t + 4), 0>,
<25, 0, [t, -t, 1], t^4*(t^2 + 4*t - 16), 0>,
<27, 1, [t, t^2, -1], t^8*(t + 2)*(t^3 - 6*t^2 - 4), 0>,
<32, 1, [t + 1, -t^3, t^2 + t], t^20*(t + 1)^6*(4*t^4 - 12*t^2 - 16*t - 7), 0 >,
<36, 1, [t + 1, -t^2 - t - 1, -t^2 - t],
  t^12*(t + 1)^4*(4*t^4 + 8*t^3 + 12*t^2 + 8*t + 1), 0>,
<49, 1, [t, -t^2, t - 1],
  t^8*(t - 2)^2*(t - 1)^4*(t + 1)^2*(t^4 - 6*t^3 + 3*t^2 + 18*t - 19), 0>
];

function SetPrecision(m)
  prec := 5;
  if m in [2, 3, 4, 5] then prec := 15; end if;
  if m in [12, 15, 16] then prec := 25; end if;
  if m in [8, 11, 14] then prec := 30; end if;
  if m in [7, 9] then prec := 50; end if;
  return prec;
end function;

GenusZeroSubsts := func<m,t| case<m|
    2 : -64*t/(t + 64)^2,
    3 : -2*t/(t + 27)^2,
    4 : -2*t/(t + 16)^2,
    5 : t/(t^2 + 22*t + 125),
    6 : -(t + 8)*(t + 9)/t,
    7 : (t^2 + 13*t + 49)/t,
    8 : -(t^2 + 14*t + 32)/t,
    9 : (t^2 + 9*t + 27)/t,
    10 : -(t^2 + 8*t + 20)/t,
    12 : -t/((t + 3)*(t + 4)),
    16 : (t^2 + 5*t + 8)/(t + 2)/(t + 4),
    18 : -t/((t + 2)*(t + 3)),
    25 : -(t^2 + 2*t + 5)/t,
    default : 0 >>;

GenusOneSubsts := func<m,u,v| case<m|
    11 : (v + 3*u + 13)/(u - 16),
    14 : -(v - 3*u + 4)/((u - 1)*(u - 2)),
    15 : (v - (u + 1)*(u + 4))/((u + 1)*(u + 2)),
    17 : (v + 2*u + 7)/(u - 7),
    19 : -(v + 2*u)/(u - 5),
    20 : -(v + (u + 1)*(u - 2))/(u*(u + 1)),
    21 : -(v + 2*u + 3)/(u - 5),
    24 : (v + u + 2)/(u - 4),
    27 : (v + u + 2)/(u - 3),
    32 : (v - u^2)/(u*(u - 2)),
    36 : (v - u^2 + 1)/(u*(u - 2)),
    49 : (v + 2*u - 3)/(u - 2),
    default : 0 >>;

print "Checking Table 1 (including the moduli interpretations)";

mm := [m : m in [2..1000] | GenusX0N(m) in [0,1] and GCD(m,13) eq 1];

for ctr in [1..#Table1] do
  m,g,subst,c,r := Explode(Table1[ctr]);
  print m,g,[Sprint(u): u in subst];
  assert m eq mm[ctr];
  assert g eq GenusX0N(m);
  assert g in [0,1];
  k := LegendreSymbol(m,13) eq 1 select 1 else 2;
  LHS := Evaluate(F[k],subst) + O(e^(r+1));
  RHS := c*e^r + O(e^(r+1));
  assert LHS eq RHS;
  b := SetPrecision(m);
  sigma := Evaluate(jsum[k],[xi + O(e^b): xi in subst]);
  tau := Evaluate(jprod[k],[xi + O(e^b): xi in subst]);
  sigma := Evaluate(sigma,0);
  tau := Evaluate(tau,0);
  X0m := SmallModularCurve(m);
  if g eq 0 then
    KK := Kt;
    X0m_genpt := BaseChange(X0m,KK)![t,1];
    t1 := GenusZeroSubsts(m,t);
  elif g eq 1 then
    KK<u,v> := FunctionField(X0m);
    X0m_genpt := BaseChange(X0m,KK)![u,v];
    t1 := GenusOneSubsts(m,u,v);
  end if;
  j := jNInvariant(X0m_genpt,m);
  j1 := jInvariant(X0m_genpt,m);
  assert IsSquare(Evaluate(c,t1));
  P<X> := PolynomialRing(KK);
  jpoly := X^2 - Evaluate(sigma,t1)*X + Evaluate(tau,t1);
  assert {r[1] : r in Roots(jpoly)} eq {j,j1};
  if m eq 10 then // example in the text
    assert j1 eq Evaluate(j,20/t);
    assert t1 eq -(t^2 + 8*t + 20)/t;
    assert j eq (t*(t + 4)^5 + 16*(t + 5))^3/(t*(t + 4)^5*(t + 5)^2);
  end if; 
end for;

print "\nThe locus j = j'";

// The formulae for the j-maps from Theorem 2.4

A1,B1,D1 := Z13ModuliInterpretation(1,x/z,y/z);
A1 := R!(z^24*A1);
B1 := R!(z^36*B1);
D1 := R!(z^72*D1);
Delta1 := x^4 + x^3*y - x^3*z + x*y^2*z - x*y*z^2 - y^2*z^2 + y*z^3;
mu1 := y^5*z^7*(x + y - z)^4*(x^2 + y*z - z^2)^2;
assert D1 eq mu1*Delta1^13;

A2,B2,D2 := Z13ModuliInterpretation(2,x/z,y/z);
A2 := R!(z^24*A2);
B2 := R!(z^36*B2);
D2 := R!(z^72*D2);
Delta2 := x^3*y + x^2*y^2 + 2*x^2*y*z + x*y^2*z + x*y*z^2 + x*z^3 + y*z^3;
mu2 := -x^6*z^8*(x^2 + x*y + x*z + z^2)^3;
assert D2 eq mu2*Delta2^13;

A := [A1,A2];
B := [B1,B2];
D := [D1,D2];
Delta := [Delta1,Delta2];
mu := [mu1,mu2];

// Some irreducible factors of the numerator of (j - j')^2

G := [ x^7*z - 4*x^6*z^2 - x^5*y^3 + x^5*y*z^2 + 6*x^5*z^3 - x^4*y^4 +
    7*x^4*y^3*z - 4*x^4*y^2*z^2 + x^4*y*z^3 - 4*x^4*z^4 + 5*x^3*y^4*z
    - 23*x^3*y^3*z^2 + 26*x^3*y^2*z^3 - 9*x^3*y*z^4 + x^3*z^5 -
    2*x^2*y^5*z - 4*x^2*y^4*z^2 + 23*x^2*y^3*z^3 - 28*x^2*y^2*z^4 +
    11*x^2*y*z^5 + 6*x*y^5*z^2 - 26*x*y^4*z^3 + 30*x*y^3*z^4 -
    6*x*y^2*z^5 - 4*x*y*z^6 - y^6*z^2 - 10*y^5*z^3 + 35*y^4*z^4 -
    36*y^3*z^5 + 12*y^2*z^6, x^11 + 3*x^10*y + 3*x^9*y^2 + 7*x^9*y*z -
    2*x^9*z^2 + x^8*y^3 + 18*x^8*y^2*z - 11*x^8*y*z^2 + 15*x^7*y^3*z +
    x^7*y^2*z^2 - 17*x^7*y*z^3 + x^7*z^4 + 4*x^6*y^4*z +
    31*x^6*y^3*z^2 - 59*x^6*y^2*z^3 + 24*x^6*y*z^4 + 27*x^5*y^4*z^2 -
    43*x^5*y^3*z^3 + 17*x^5*y^2*z^4 - 6*x^5*y*z^5 + 6*x^4*y^5*z^2 +
    17*x^4*y^4*z^3 - 57*x^4*y^3*z^4 + 35*x^4*y^2*z^5 + 21*x^3*y^5*z^3
    - 54*x^3*y^4*z^4 + 51*x^3*y^3*z^5 - 18*x^3*y^2*z^6 + 4*x^2*y^6*z^3
    - 2*x^2*y^5*z^4 - 7*x^2*y^4*z^5 + 5*x^2*y^3*z^6 + 6*x*y^6*z^4 -
    16*x*y^5*z^5 + 14*x*y^4*z^6 - 4*x*y^3*z^7 + y^7*z^4 - 3*y^6*z^5 +
    3*y^5*z^6 - y^4*z^7, x^6*y*z^2 + x^5*y^4 + 11*x^5*y*z^3 +
    2*x^4*y^5 + 6*x^4*y^4*z + 2*x^4*y^2*z^3 + 35*x^4*y*z^4 + x^4*z^5 +
    x^3*y^6 + 6*x^3*y^5*z + 8*x^3*y^4*z^2 + 2*x^3*y^3*z^3 +
    20*x^3*y^2*z^4 + 37*x^3*y*z^5 + 10*x^3*z^6 + 3*x^2*y^5*z^2 +
    15*x^2*y^4*z^3 + 10*x^2*y^3*z^4 + 2*x^2*y^2*z^5 + 41*x^2*y*z^6 +
    10*x^2*z^7 + 3*x*y^4*z^4 + 12*x*y^3*z^5 + 7*x*y^2*z^6 + 8*x*y*z^7
    + 10*x*z^8 + y^3*z^6 + 3*y^2*z^7 + 3*y*z^8 + z^9, x^8*y^2 +
    2*x^7*y^3 + 10*x^7*y^2*z + x^6*y^4 + 10*x^6*y^3*z + 39*x^6*y^2*z^2
    + 2*x^6*y*z^3 + 31*x^5*y^3*z^2 + 63*x^5*y^2*z^3 + 18*x^5*y*z^4 +
    29*x^4*y^3*z^3 + 56*x^4*y^2*z^4 + 48*x^4*y*z^5 + x^4*z^6 +
    36*x^3*y^2*z^5 + 49*x^3*y*z^6 + 8*x^3*z^7 - 6*x^2*y^2*z^6 +
    35*x^2*y*z^7 + 10*x^2*z^8 + 3*x*y*z^8 + 8*x*z^9 + y*z^9 + z^10,
    x^9*y^2 + 4*x^8*y^3 + 8*x^8*y^2*z + 6*x^7*y^4 + 24*x^7*y^3*z +
    35*x^7*y^2*z^2 + 2*x^7*y*z^3 + 4*x^6*y^5 + 24*x^6*y^4*z +
    71*x^6*y^3*z^2 + 87*x^6*y^2*z^3 + 14*x^6*y*z^4 + x^5*y^6 +
    8*x^5*y^5*z + 37*x^5*y^4*z^2 + 112*x^5*y^3*z^3 + 140*x^5*y^2*z^4 +
    46*x^5*y*z^5 + x^5*z^6 + x^4*y^5*z^2 + 27*x^4*y^4*z^3 +
    98*x^4*y^3*z^4 + 143*x^4*y^2*z^5 + 82*x^4*y*z^6 + 6*x^4*z^7 +
    35*x^3*y^3*z^5 + 92*x^3*y^2*z^6 + 82*x^3*y*z^7 + 12*x^3*z^8 +
    21*x^2*y^2*z^7 + 42*x^2*y*z^8 + 15*x^2*z^9 + 7*x*y*z^9 + 7*x*z^10
    + z^11 ];
    
H := [ x^8 + 2*x^7*y + 2*x^7*z + x^6*y^2 + x^6*y*z - x^6*z^2 +
    2*x^5*y^2*z + 12*x^5*y*z^2 - 4*x^5*z^3 + 3*x^4*y^3*z +
    9*x^4*y^2*z^2 - 6*x^4*y*z^3 - x^4*z^4 - 2*x^3*y^3*z^2 +
    19*x^3*y^2*z^3 - 25*x^3*y*z^4 + 2*x^3*z^5 + 3*x^2*y^4*z^2 +
    15*x^2*y^3*z^3 - 32*x^2*y^2*z^4 + 14*x^2*y*z^5 + x^2*z^6 -
    2*x*y^4*z^3 + 10*x*y^3*z^4 - 19*x*y^2*z^5 + 11*x*y*z^6 + y^5*z^3 +
    6*y^4*z^4 - 23*y^3*z^5 + 25*y^2*z^6 - 9*y*z^7, x^11 + 5*x^10*y -
    x^10*z + 10*x^9*y^2 - 3*x^9*y*z - 2*x^9*z^2 + 10*x^8*y^3 +
    x^8*y^2*z - 20*x^8*y*z^2 + 2*x^8*z^3 + 5*x^7*y^4 + 14*x^7*y^3*z -
    45*x^7*y^2*z^2 + 21*x^7*y*z^3 + x^7*z^4 + x^6*y^5 + 21*x^6*y^4*z -
    33*x^6*y^3*z^2 + 2*x^6*y^2*z^3 + 11*x^6*y*z^4 - x^6*z^5 +
    13*x^5*y^5*z + 4*x^5*y^4*z^2 - 66*x^5*y^3*z^3 + 69*x^5*y^2*z^4 -
    20*x^5*y*z^5 + 3*x^4*y^6*z + 18*x^4*y^5*z^2 - 59*x^4*y^4*z^3 +
    61*x^4*y^3*z^4 - 29*x^4*y^2*z^5 + 6*x^4*y*z^6 + 11*x^3*y^6*z^2 -
    4*x^3*y^5*z^3 - 40*x^3*y^4*z^4 + 48*x^3*y^3*z^5 - 15*x^3*y^2*z^6 +
    3*x^2*y^7*z^2 + 8*x^2*y^6*z^3 - 46*x^2*y^5*z^4 + 63*x^2*y^4*z^5 -
    35*x^2*y^3*z^6 + 7*x^2*y^2*z^7 + 3*x*y^7*z^3 - 3*x*y^6*z^4 -
    9*x*y^5*z^5 + 15*x*y^4*z^6 - 6*x*y^3*z^7 + y^8*z^3 + y^7*z^4 -
    13*y^6*z^5 + 23*y^5*z^6 - 16*y^4*z^7 + 4*y^3*z^8, x^13 + 3*x^12*y
    - 5*x^12*z + 3*x^11*y^2 - 5*x^11*y*z + 10*x^11*z^2 + x^10*y^3 +
    10*x^10*y^2*z - 8*x^10*y*z^2 - 10*x^10*z^3 + 15*x^9*y^3*z -
    25*x^9*y^2*z^2 + 22*x^9*y*z^3 + 5*x^9*z^4 + 5*x^8*y^4*z +
    13*x^8*y^3*z^2 - 10*x^8*y^2*z^3 - 13*x^8*y*z^4 - x^8*z^5 +
    30*x^7*y^4*z^2 - 55*x^7*y^3*z^3 + 33*x^7*y^2*z^4 - x^7*y*z^5 +
    10*x^6*y^5*z^2 + 8*x^6*y^4*z^3 - 13*x^6*y^3*z^4 - 2*x^6*y^2*z^5 +
    2*x^6*y*z^6 + 30*x^5*y^5*z^3 - 71*x^5*y^4*z^4 + 46*x^5*y^3*z^5 -
    11*x^5*y^2*z^6 + 10*x^4*y^6*z^3 + x^4*y^5*z^4 - 15*x^4*y^4*z^5 +
    3*x^4*y^3*z^6 + 2*x^4*y^2*z^7 + 15*x^3*y^6*z^4 - 46*x^3*y^5*z^5 +
    39*x^3*y^4*z^6 - 8*x^3*y^3*z^7 + 5*x^2*y^7*z^4 - 2*x^2*y^6*z^5 -
    11*x^2*y^5*z^6 + 10*x^2*y^4*z^7 - 2*x^2*y^3*z^8 + 3*x*y^7*z^5 -
    11*x*y^6*z^6 + 13*x*y^5*z^7 - 5*x*y^4*z^8 + y^8*z^5 - y^7*z^6 -
    2*y^6*z^7 + 3*y^5*z^8 - y^4*z^9, x^6*y + 3*x^5*y^2 + 3*x^4*y^3 +
    2*x^4*y*z^2 + x^4*z^3 + x^3*y^4 - 4*x^3*y^2*z^2 + 13*x^3*y*z^3 -
    x^3*z^4 - 6*x^2*y^3*z^2 + 5*x^2*y^2*z^3 + 7*x^2*y*z^4 + x*y^3*z^3
    - 4*x*y^2*z^4 - x*y*z^5 + 8*x*z^6 + y^2*z^5 + 2*y*z^6 + z^7, x^9*y
    + 4*x^8*y^2 + 6*x^7*y^3 - 4*x^7*y*z^2 + x^7*z^3 + 4*x^6*y^4 -
    8*x^6*y^2*z^2 - x^6*y*z^3 - x^6*z^4 + x^5*y^5 - 4*x^5*y^3*z^2 +
    6*x^5*y*z^4 - 4*x^5*z^5 + 9*x^4*y^3*z^3 + 10*x^4*y^2*z^4 +
    19*x^4*y*z^5 - 2*x^4*z^6 + 7*x^3*y^4*z^3 + 7*x^3*y^3*z^4 +
    35*x^3*y^2*z^5 + 9*x^3*y*z^6 + 11*x^3*z^7 + 23*x^2*y^3*z^5 +
    9*x^2*y^2*z^6 + 12*x^2*y*z^7 + 14*x^2*z^8 + x*y^3*z^6 +
    19*x*y^2*z^7 - 7*x*y*z^8 + 11*x*z^9 + y^2*z^8 + 2*y*z^9 + z^10,
    x^10*y^2 + 3*x^9*y^3 + 4*x^9*y^2*z + 3*x^8*y^4 + 8*x^8*y^3*z +
    17*x^8*y^2*z^2 + 2*x^8*y*z^3 + x^7*y^5 + 4*x^7*y^4*z +
    32*x^7*y^3*z^2 + 45*x^7*y^2*z^3 + 6*x^7*y*z^4 + 15*x^6*y^4*z^2 +
    61*x^6*y^3*z^3 + 77*x^6*y^2*z^4 + 22*x^6*y*z^5 + x^6*z^6 +
    18*x^5*y^4*z^3 + 65*x^5*y^3*z^4 + 113*x^5*y^2*z^5 + 44*x^5*y*z^6 +
    2*x^5*z^7 + 58*x^4*y^3*z^5 + 97*x^4*y^2*z^6 + 65*x^4*y*z^7 +
    6*x^4*z^8 + x^3*y^3*z^6 + 69*x^3*y^2*z^7 + 57*x^3*y*z^8 +
    10*x^3*z^9 + 3*x^2*y^2*z^8 + 36*x^2*y*z^9 + 11*x^2*z^10 +
    3*x*y*z^10 + 7*x*z^11 + z^12 ];

M := x^8 + 4*x^7*y - 4*x^7*z + 6*x^6*y^2 - 21*x^6*y*z + 6*x^6*z^2 +
    4*x^5*y^3 - 37*x^5*y^2*z + 43*x^5*y*z^2 - 4*x^5*z^3 + x^4*y^4 -
    25*x^4*y^3*z + 68*x^4*y^2*z^2 - 43*x^4*y*z^3 + x^4*z^4 -
    3*x^3*y^4*z + 22*x^3*y^3*z^2 - 28*x^3*y^2*z^3 + 21*x^3*y*z^4 +
    2*x^2*y^5*z - 17*x^2*y^4*z^2 + 63*x^2*y^3*z^3 - 39*x^2*y^2*z^4 -
    4*x^2*y*z^5 - 7*x*y^5*z^2 + 65*x*y^4*z^3 - 103*x*y^3*z^4 +
    39*x*y^2*z^5 + y^6*z^2 + 13*y^5*z^3 - 43*y^4*z^4 + 39*y^3*z^5 -
    9*y^2*z^6;

assert [Degree(f): f in G] eq [8,11,9,10,11];
assert [Degree(f): f in H] eq [8,11,13,7,10,12];

// Checking that the numerators and denominators of (j - j')^2
// factor as claimed in the paper

k := 1;
jprod := A[k]^3/D[k];
jsum := (jprod + 1728^2 - B[k]^2/D[k])/1728;
num := F[1]*G[1]^2*G[2]^2*H[1]^2*H[2]^2*H[3]^2*M^2; 
den := y^8*z^8*(x + y - z)^6*(x^2 + y*z - z^2)^2*Delta[1]^26;
assert jsum^2 - 4*jprod eq num/den;

k := 2;
jprod := A[k]^3/D[k];
jsum := (jprod + 1728^2 - B[k]^2/D[k])/1728;
num := F[2]*G[3]^2*G[4]^2*G[5]^2*H[4]^2*H[5]^2*H[6]^2;
den := x^10*z^10*(x^2 + x*y + x*z + z^2)^4*Delta[2]^26;
assert jsum^2 - 4*jprod eq num/den;

// Baran's modular curve of level 13

Beqn := (y + z)*x^3 - (2*y^2 + y*z)*x^2
         + (y^3 - y^2*z + 2*y*z^2 - z^3)*x - (2*y^2*z^2 - 3*y*z^3);
X13B := Curve(P2,Beqn);

// Formulae for the j-map matching Appendix A of Baran's paper.
// (In fact computed as described in Section 2.3)

function jMapBaran(label,X,Y)
  assert label in ["s","ns"];
  if label eq "s" then 
    g_s := (-Y^44 - 478*Y^43 - 32789*Y^42 - 762363*Y^41 - 7547176*Y^40
      - 29112957*Y^39 + 27142557*Y^38 + 516914701*Y^37 +
      812780834*Y^36 - 3217396462*Y^35 - 9372779891*Y^34 +
      10412400391*Y^33 + 51281877696*Y^32 - 21282455533*Y^31 -
      192026183004*Y^30 + 36250779915*Y^29 + 569052119000*Y^28 -
      70791418862*Y^27 - 1392885057099*Y^26 + 180780916888*Y^25 +
      2828775392943*Y^24 - 517934371345*Y^23 - 4803267109749*Y^22 +
      1263736631868*Y^21 + 6827480292072*Y^20 - 2364214517403*Y^19 -
      8006349554252*Y^18 + 3419622432067*Y^17 + 7623483793639*Y^16 -
      3825699399717*Y^15 - 5822086727508*Y^14 + 3183399803606*Y^13 +
      3469776050052*Y^12 - 1874600063602*Y^11 - 1524399772882*Y^10 +
      753079641994*Y^9 + 455793324718*Y^8 - 201254434444*Y^7 -
      83894215212*Y^6 + 34538107094*Y^5 + 8027212027*Y^4 -
      3468719981*Y^3 - 204764859*Y^2 + 159120623*Y - 12145116)*X^2 +
      (Y^45 + 446*Y^44 + 27904*Y^43 + 571184*Y^42 + 4618918*Y^41 +
      10200852*Y^40 - 54213498*Y^39 - 267763439*Y^38 + 165832967*Y^37
      + 2368932904*Y^36 + 395119701*Y^35 - 12956472125*Y^34 -
      4150166662*Y^33 + 53633145847*Y^32 + 12575521646*Y^31 -
      180183851743*Y^30 - 10692143892*Y^29 + 494822838826*Y^28 -
      75462285877*Y^27 - 1111012707457*Y^26 + 427597488033*Y^25 +
      2035065625744*Y^24 - 1283704039956*Y^23 - 2989648304276*Y^22 +
      2716215586599*Y^21 + 3407000739960*Y^20 - 4339306349468*Y^19 -
      2827171518873*Y^18 + 5300271455330*Y^17 + 1407694779775*Y^16 -
      4900196317351*Y^15 + 18526414704*Y^14 + 3345715118980*Y^13 -
      643829576586*Y^12 - 1623126310814*Y^11 + 509263027702*Y^10 +
      527959900322*Y^9 - 203278262234*Y^8 - 104991650269*Y^7 +
      45780617674*Y^6 + 10651697173*Y^5 - 5549670458*Y^4 -
      270006722*Y^3 + 269295673*Y^2 - 26210553*Y - 144612)*X - 2*Y^44
      - 891*Y^43 - 55361*Y^42 - 1113897*Y^41 - 8622455*Y^40 -
      14660030*Y^39 + 130330785*Y^38 + 526492528*Y^37 - 657396776*Y^36
      - 5427279992*Y^35 + 569333662*Y^34 + 32403186745*Y^33 +
      8217537249*Y^32 - 137416119935*Y^31 - 44405672291*Y^30 +
      459234305539*Y^29 + 116561534270*Y^28 - 1261675920964*Y^27 -
      154320890766*Y^26 + 2870046210200*Y^25 - 102772109156*Y^24 -
      5399242483910*Y^23 + 1112801469826*Y^22 + 8364611871299*Y^21 -
      3167545239322*Y^20 - 10526258387861*Y^19 + 5800415621050*Y^18 +
      10505818581326*Y^17 - 7693864828631*Y^16 - 8045599613602*Y^15 +
      7523007777896*Y^14 + 4507918282558*Y^13 - 5289650276066*Y^12 -
      1710818708766*Y^11 + 2566391186088*Y^10 + 370366462290*Y^9 -
      818009316691*Y^8 - 14872094943*Y^7 + 160086948432*Y^6 -
      12957280360*Y^5 - 16931385812*Y^4 + 2913054190*Y^3 +
      664153911*Y^2 - 193921364*Y + 12288000;
    d_s := X^2*Y + X^2 - X - 1;
    j := g_s/d_s^13;
  else
    g_ns := (884736*Y^55 + 49489920*Y^54 + 1012211712*Y^53 +
      11898450496*Y^52 + 95148841568*Y^51 + 564433986856*Y^50 +
      2616961812219*Y^49 + 9822086693987*Y^48 + 30642386133997*Y^47 +
      81359505328755*Y^46 + 188611532400778*Y^45 +
      392863091096565*Y^44 + 750115492349610*Y^43 +
      1285021297306328*Y^42 + 1751719154883360*Y^41 +
      1150764689309236*Y^40 - 2039149007306334*Y^39 -
      5717108163092124*Y^38 + 7099555968682235*Y^37 +
      77776536294315788*Y^36 + 240934894546612597*Y^35 +
      395125111510297175*Y^34 + 89165142545864768*Y^33 -
      1550321042833254477*Y^32 - 5305245647359351740*Y^31 -
      10446908162336304678*Y^30 - 12871375409731491993*Y^29 -
      4395668686472870357*Y^28 + 24906193210193909056*Y^27 +
      81237769117029166717*Y^26 + 161268611204135399411*Y^25 +
      249955252927207983657*Y^24 + 324483552163739360618*Y^23 +
      363683177940575393974*Y^22 + 358040773388507611446*Y^21 +
      313163052354472751348*Y^20 + 244627586482873074382*Y^19 +
      170321443968359565942*Y^18 + 105533439538375007957*Y^17 +
      59274776728757700369*Y^16 + 31575303637834908430*Y^15 +
      16108641614130173307*Y^14 + 6982396840373322235*Y^13 +
      1946153985121206304*Y^12 + 125999715600544888*Y^11 -
      17233649727426000*Y^10 + 97243737410324752*Y^9 +
      60793442160963072*Y^8 + 4429798119424768*Y^7 -
      5098133689396992*Y^6 - 508472036045568*Y^5 + 376022069526528*Y^4
      + 11471461152768*Y^3 - 15545480048640*Y^2 + 1852894310400*Y -
      67903488000)*X^2 + (-884736*Y^56 - 41306112*Y^55 -
      745311744*Y^54 - 7854617152*Y^53 - 56752908832*Y^52 -
      305691731184*Y^51 - 1291334290291*Y^50 - 4426288172137*Y^49 -
      12647728750928*Y^48 - 31075240957232*Y^47 - 69032547442536*Y^46
      - 148896176655469*Y^45 - 322250360417387*Y^44 -
      632799187856807*Y^43 - 813902597999889*Y^42 +
      446780681983306*Y^41 + 5558180767681622*Y^40 +
      12861196240238133*Y^39 + 82777469121694*Y^38 -
      95303288589756575*Y^37 - 346077988812264244*Y^36 -
      660107106734655623*Y^35 - 436018183788282799*Y^34 +
      1681210077303367073*Y^33 + 7273388110177545277*Y^32 +
      16144187701672166507*Y^31 + 23113156489251424786*Y^30 +
      15841883577593033385*Y^29 - 22985886157655295827*Y^28 -
      107891300530866582629*Y^27 - 240236673343003218988*Y^26 -
      401070259854975037146*Y^25 - 553221859660242994706*Y^24 -
      655402208407230091461*Y^23 - 681436550533356290765*Y^22 -
      630169426685781559948*Y^21 - 521000422369504316094*Y^20 -
      384796397343162994888*Y^19 - 255121295576640474717*Y^18 -
      155328774650536353143*Y^17 - 88765501587983950847*Y^16 -
      45750129339254668296*Y^15 - 18659682247552136243*Y^14 -
      5068590639242514778*Y^13 - 1175381774283273326*Y^12 -
      905896315537604416*Y^11 - 627696074121250176*Y^10 -
      142339763530459808*Y^9 + 40822149606469600*Y^8 +
      17042318686921728*Y^7 - 4711676288999168*Y^6 -
      1772235425484288*Y^5 + 359764019596800*Y^4 + 59190591062016*Y^3
      - 21787540807680*Y^2 + 2089099468800*Y - 67931136000)*X +
      1769472*Y^55 + 81727488*Y^54 + 1448432640*Y^53 +
      14896074368*Y^52 + 104157512704*Y^51 + 536681652096*Y^50 +
      2133314691734*Y^49 + 6717439350959*Y^48 + 17014862466548*Y^47 +
      35253668331292*Y^46 + 62993289306098*Y^45 + 112164904174479*Y^44
      + 225104747950936*Y^43 + 392976024258106*Y^42 +
      41136059444059*Y^41 - 3005798979372973*Y^40 -
      11839374362691297*Y^39 - 21190199356139685*Y^38 +
      6071591688951686*Y^37 + 156325146877196166*Y^36 +
      498229160275748461*Y^35 + 823122591905419376*Y^34 +
      242278135736846249*Y^33 - 2876324971563594929*Y^32 -
      9820974562808688564*Y^31 - 19019630528554710562*Y^30 -
      23294749281381561713*Y^29 - 9486631575477114590*Y^28 +
      36945136929714190121*Y^27 + 122510790902234508091*Y^26 +
      237163604354617256618*Y^25 + 354702102863369258879*Y^24 +
      443147794960495122187*Y^23 + 477688281307028925796*Y^22 +
      450443978678593045613*Y^21 + 374770111259894186150*Y^20 +
      278602970647572826148*Y^19 + 187430241355450245572*Y^18 +
      113084138555091016698*Y^17 + 58445464663933323835*Y^16 +
      24688966627452318352*Y^15 + 9359950076743879942*Y^14 +
      4515705377792888754*Y^13 + 2518265993855063344*Y^12 +
      865069239507223888*Y^11 - 38745965042080*Y^10 -
      90275758453065632*Y^9 + 96762513188096*Y^8 +
      14322780344679936*Y^7 + 1383116431053312*Y^6 -
      986975640101376*Y^5 - 23199406829568*Y^4 + 43522275594240*Y^3 -
      5384442470400*Y^2 + 203793408000*Y;
    d_ns := -(Y + 2)*X^2 + (Y^2 + 5*Y)*X - (Y^2 - Y);
    j := g_ns/((Y + 1)^24*d_ns^13);
  end if;
  return j;
end function;

// The double covers X_s(13) -> X^+_s(13) and X_ns(13) -> X^+_ns(13)
// each correspond to a degree 2 extension of function fields
// with Kummer generators delta_s and delta_ns as follows

KB<u,v> := FunctionField(X13B);
assert Evaluate(Beqn,[u,v,1]) eq 0;

delta_s := 5*u^2*v^2 - 2*u^2*v - 7*u^2 - 4*u*v^3 + 10*u*v^2 -
  4*u*v - 6*u + 5*v^2 - 10*v - 3;

delta_ns := -8*u^4 + 20*u^3*v + 20*u^2*v^2 - 44*u^2*v + u^2 -
  40*u*v^3 + 114*u*v^2 - 52*u*v - 11*v^4 - 64*v^3 + 80*v^2;

// Formulae to help identify the curves defined by the F_i,G_i,H_i.

substs := [ [ 20*u^4 - 99*u^3*v + 7*u^3 + 43*u^2*v^2 + 186*u^2*v -
    2*u^2 - 14*u*v^3 - 109*u*v^2 + 162*u*v - 124*u - 8*v^4 + 28*v^3 +
    47*v^2 - 30*v - 126, 10*u^4 - 26*u^3*v + 117*u^3 + 101*u^2*v^2 -
    427*u^2*v + 38*u^2 - 19*u*v^3 + 173*u*v^2 - 189*u*v + 183*u -
    8*v^4 + 76*v^3 - 215*v^2 + 3*v + 252, -30*u^4 - 34*u^3*v + 17*u^3
    + 142*u^2*v^2 - 85*u^2*v + 51*u^2 - 62*u*v^3 + 46*u*v^2 + 10*u*v -
    50*u - 16*v^4 + 104*v^3 - 110*v^2 + 12*v - 63 ], [ 12*u^4 +
    51*u^3*v^2 + 14*u^3*v + 7*u^3 - 93*u^2*v^3 - 28*u^2*v^2 + 24*u^2*v
    - 5*u^2 + 45*u*v^4 - 37*u*v^3 - 4*u*v^2 - 2*u - 8*v^5 - 9*v^4 -
    23*v^3 + 101*v^2 + 14*v, -9*u^4 - 40*u^3*v^2 - 14*u^3*v - 7*u^3 +
    66*u^2*v^3 + 7*u^2*v^2 - 18*u^2*v + 2*u^2 - 16*u*v^4 + 63*u*v^3 -
    3*u*v^2 + 6*u*v + 2*u + 32*v^3 - 78*v^2 - 11*v, -6*u^4 -
    39*u^3*v^2 - 12*u^3*v - 7*u^3 + 71*u^2*v^3 - 7*u^2*v^2 - 10*u^2*v
    - u^2 - 27*u*v^4 + 71*u*v^3 - 9*u*v^2 + 10*u*v + 2*u - 16*v^4 +
    20*v^3 - 69*v^2 - 10*v ], [ 4*u^2*v + 2*u^2 - 2*u*v^2 + v^3 + v^2
    - 3*v - 2, 2*u^2*v + 2*u^2 - u*v^2 + 2*u*v - v^3 + 3*v - 2, -u^3 +
    2*u^2*v + u^2 - u*v^2 + u*v + 4*u - v^2 - v + 2 ], [ u^3 -
    2*u^2*v^2 - 7*u^2*v + 4*u*v^3 + 2*u*v^2 - 5*u*v - u - 2*v^4 +
    9*v^3 - 4*v^2 - 3*v, -4*u^3*v - 7*u^3 + u^2*v^2 + 12*u^2*v - 6*u^2
    + 6*u*v^3 - 12*u*v^2 + 11*u*v + 7*u - 8*v^2 - 4*v + 6, -3*u^3*v -
    2*u^3 + u^2*v^2 - 4*u^2*v + 3*u*v^3 - 5*u*v + 6*u + 2*v^4 - v^3 -
    4*v^2 - v + 6 ], [ 2*u^3 - u^2*v - u*v^2 + 3*u*v - 2*u + v^3 -
    2*v^2 + 2*v, -u^2*v - u^2 + 2*u*v^2 - u*v - 2*u - v^3 + 2*v^2 -
    2*v, -u^2*v + u^2 + 2*u*v^2 - 2*u*v + 2*u - v^3 + 2*v^2 - 2*v ], [
    u^4 + u^3*v - 2*u^3 + u^2*v^2 + u^2*v - u*v^3 - 4*u*v + 2*u + v^2
    + 2*v - 1, -u^4 - 2*u^3*v - u^3 - u^2*v^2 - 9*u^2*v + 2*u*v^3 +
    2*u*v^2 - 4*u*v + u - 3*v^2 + v + 1, -2*u^4 - 2*u^3*v - u^3 +
    2*u^2*v^2 - 3*u^2*v + 4*u^2 + u*v^2 - u*v + u - v^2 - 2 ], [ u^3 -
    3*u^2*v + u^2 + 3*u*v^2 - 3*u*v + u - v^3 + 2*v^2 - v + 1, u^3*v -
    u^3 - 3*u^2*v^2 + 3*u^2*v - u^2 + 2*u*v - v^3 + 3*v^2 - 2*v - 1,
    -u^3*v + 2*u^2*v^2 - u*v^3 + 2*u*v - u + v^3 - 2*v^2 + v - 1 ], [
    u^3 + 3*u^2 + 2*u - v^2, u^3 - 4*u^2*v + u^2 + u*v^2 - 3*u*v + u +
    v^3 - 2*v^2 - 3*v, 2*u^3 - 4*u^2*v + 4*u^2 - 3*u*v + 4*u + v^3 -
    2*v^2 + 2 ], [ 2*u^3*v - 4*u^3 + 6*u^2*v^2 + 43*u^2*v + 25*u^2 -
    7*u*v^3 - 27*u*v^2 - 14*u*v + 4*u - 13*v^4 - 32*v^3 + 13*v^2 +
    36*v, -3*u^4 + 11*u^3*v + 3*u^3 - 13*u^2*v^2 - 21*u^2*v + 3*u^2 -
    20*u*v^3 - 33*u*v - 3*u + 13*v^4 - 7*v^3 - 16*v^2 - 7*v, -3*u^4 +
    4*u^3*v - 4*u^3 - 13*u^2*v^2 - 16*u^2*v + 3*u^2 + 13*u*v^3 +
    26*u*v^2 + 15*u*v + 29*u - 13*v^4 - 45*v^3 - 59*v^2 - 51*v ], [
    19*u^3*v + 293*u^3 + 297*u^2*v^2 - 666*u^2*v + 360*u^2 - 540*u*v^3
    + 515*u*v^2 - 589*u*v - 113*u + 224*v^4 - 389*v^3 + 74*v^2 + 629*v
    + 41, -192*u^3*v - 192*u^3 + 608*u^2*v^2 + 588*u^2*v - 376*u^2 -
    640*u*v^3 - 824*u*v^2 + 551*u*v - 389*u + 224*v^4 + 507*v^3 +
    62*v^2 + 209*v - 205, 137*u^4 - 570*u^3*v + 375*u^3 + 1062*u^2*v^2
    - 713*u^2*v + 223*u^2 - 966*u*v^3 + 106*u*v^2 - 61*u*v - 195*u +
    448*v^4 + 118*v^3 - 292*v^2 + 670*v + 41 ], [ u*v^2 + u*v - v,
    u^2*v + 3*u^2 - 2*u*v^2 - 6*u*v + 3*u + v^3 + 2*v^2 - 4*v, u*v^2 +
    u*v - v^3 + 2*v ], [ 4*u^4 - 9*u^3*v + 5*u^3 + 5*u^2*v^2 - 6*u^2*v
    + 2*u^2 + u*v^3 - 4*u*v^2 - u*v + u - v^4 + 4*v^3 - 2*v^2 - 2*v,
    -u^3*v - u^3 - 2*u^2*v^2 - 4*u^2*v - u^2 + 3*u*v^3 + 8*u*v^2 +
    4*u*v + 5*u + v^4 - 4*v^2 - 8*v, u^3*v - 3*u^3 - 3*u^2*v^2 +
    5*u^2*v - 3*u^2 + u*v^3 + 2*u*v^2 + 4*u*v - 2*v^3 + 2*v^2 - v ], [
    4*u^4 - 16*u^3*v + 11*u^3 + 23*u^2*v^2 - 43*u^2*v - 5*u^2 -
    25*u*v^3 + 36*u*v^2 - 14*u*v - 22*u + 3*v^4 - 34*v^3 + 34*v^2 +
    48*v, -7*u^3*v - 3*u^3 + 11*u^2*v^2 + u^2*v + 10*u*v^3 - 13*u*v^2
    - 2*u*v + 3*u - 3*v^4 + 34*v^3 - 13*v^2 - 7*v, -4*u^4 + 16*u^3*v +
    u^3 - 20*u^2*v^2 - 12*u^2*v + 4*u^2 + 19*u*v^3 + 19*u*v^2 - 22*u*v
    - u - 3*v^3 - 4*v^2 - 9*v ] ];
substs := [[f[i]/f[3]: i in [1..3]]: f in substs];

eqns := F cat G cat H;
names := ["F1","F2","G1","G2","G3","G4","G5",
          "H1","H2","H3","H4","H5","H6"];

// Checking the curves defined by the F_i,G_i,H_i have the
// moduli interpretations claimed in Sections 2.3 and 7.

for ctr in [1..#eqns] do
  f := eqns[ctr];
  soln := substs[ctr];
  assert Evaluate(f,soln) eq 0;
  k := ctr in [1,3,4,8,9,10] select 1 else 2;
  A0 := Evaluate(A[k],soln);
  B0 := Evaluate(B[k],soln);
  Delta0 := Evaluate(Delta[k],soln);
  D0 := Evaluate(mu[k],soln)*Delta0^13;
  jprod := A0^3/D0;
  jsum := (jprod + 1728^2 - B0^2/D0)/1728;
  j := jsum/2;
  assert jprod eq j^2;
  ss := [str : str in ["s","ns"] | j eq jMapBaran(str,u,v)];
  assert #ss eq 1;
  str := ss[1];
  delta := Evaluate(F[k],soln);
  if delta eq 0 then
    msg := "delta = 0";
  else
    flag1 := IsSquare(delta/delta_s);
    flag2 := IsSquare(delta/delta_ns);
    assert flag1 eq (str eq "s");
    assert flag2 eq (str eq "ns");
    msg := "";
  end if;
  printf "k = %o  Deg(%o) = %o  %o  %o\n",
    k,names[ctr],Degree(f),str,msg;
end for;

// The degree 8 polynomial M factors as the product of two quartics

K<b> := QuadraticField(13);
RK<x,y,z> := PolynomialRing(K,3);

M1 := x^4 + 2*x^3*y - 2*x^3*z + x^2*y^2 + 1/2*(-3*b - 13)*x^2*y*z +
   x^2*z^2 + 1/2*(-3*b - 7)*x*y^2*z + 1/2*(5*b + 13)*x*y*z^2 + y^3*z +
   1/2*(5*b + 13)*y^2*z^2 + (-b - 2)*y*z^3;
M2 := x^4 + 2*x^3*y - 2*x^3*z + x^2*y^2 + 1/2*(3*b - 13)*x^2*y*z +
   x^2*z^2 + 1/2*(3*b - 7)*x*y^2*z + 1/2*(-5*b + 13)*x*y*z^2 + y^3*z +
   1/2*(-5*b + 13)*y^2*z^2 + (b - 2)*y*z^3;
assert RK!M eq M1*M2;

KT<T> := PolynomialRing(K);
param := [
    6*T^3 + (-b - 7)*T - b - 1,
    (2*b - 4)*T^2 + (-2*b + 4)*T + b + 7,
    2*T^4 + (b + 3)*T^3 + (b + 1)*T^2 + (-b + 1)*T - b + 1
];
assert Evaluate(M1,param) eq 0;

// Checking that the inverse image in Z(13,1) is a copy of X_1(13)
// defined over Q(zeta_13)

f := -4*T^6 + 12*T^5 + 1/2*(11*b - 61)*T^4 + (-11*b + 41)*T^3 +
       1/2*(11*b - 61)*T^2 + 12*T - 4;
assert IsSquare(Evaluate(F[1],param)/f);
L<z> := CyclotomicField(13);
rt13 := 2*z^11 + 2*z^8 + 2*z^7 + 2*z^6 + 2*z^5 + 2*z^2 + 1;
Embed(K,L,rt13);
LT<T> := PolynomialRing(L);
f := LT!Eltseq(f);
// X_1(13) has equation y^2 = g(x) where
g := (T^3 - 3*T + 1)^2 - 2*(T^3 - 3*T + 1)*T*(T - 1) + 5*T^2*(T - 1)^2;
pi := z^11 + z^10 - z^9 + z^7 + z^6 - z^4 + z^3 + z^2;
la := 15*z^11 + 42*z^10 + 31*z^9 + 3*z^8 + 7*z^7 + 35*z^6
       + 39*z^5 + 11*z^4 + 27*z^2 + 42*z + 21;
T1 := (pi*T - (pi - 1))/((pi - 1)*T + 1);
assert ((pi - 1)*T + 1)^6*Evaluate(g,T1) eq la^2*f;

print "Group theoretic explanation";

GL2 := GL(2,13);
PGL2,q := quo<GL2|Centre(GL2)>;
PSL2 := sub<PGL2|q(GL2![0,1,-1,0]),q(GL2![1,1,0,1])>;
Cs := sub<PSL2|q(GL2![2,0,0,7])>;
Csplus := Normalizer(PSL2,Cs);
Cns := sub<PSL2|q(GL2![3,4,2,3])>;
Cnsplus := Normalizer(PSL2,Cns);
B := sub<PSL2|q(GL2![1,1,0,1])>;

gps := [Csplus,Cnsplus,Cs,Cns,B,PSL2];
names := ["Cs+","Cns+","Cs","Cns","B","PSL"];

gg := {};
for det in [1,-1] do
  printf "det = %o : ",det;
  mylist := [];
  for g in PGL2 do
    if g notin gg and ((g in PSL2) eq (det eq 1)) then
      orbit := {h*g*h^(-1): h in PSL2};
      stab := sub<PSL2|[h : h in PSL2 | h*g*h^(-1) eq g]>;
      assert exists(i){i : i in [1..#gps] | IsConjugate(PSL2,gps[i],stab)};
      Append(~mylist,i);
      gg join:= orbit;
    end if;
  end for;
  print [names[i] : i in Sort(mylist)];
end for;

print "\nSection 8\n";

CD := CremonaDatabase();
CR := func<E|Conductor(E) lt LargestConductor(CD) select
CremonaReference(E) else Sprint(Conductor(E)) cat "*">;
P<X> := PolynomialRing(Rationals());
aa := func<E|Vector(GF(13),TracesOfFrobenius(E,1000))>;

print "Example 8.3";

r := 4/3;
s := 5/3;
g,h := Z13Equations(1,r,s);
assert IsSquare(h^2 + 4*g);
A,B,D := Z13ModuliInterpretation(1,r,s);
jprod := A^3/D;
jsum := (jprod + 1728^2 - B^2/D)/1728;
jj := {r[1] : r in Roots(X^2 - jsum*X + jprod)};
E1 := EllipticCurve([1,0,1,-464619,-122105558]);
E2 := EllipticCurve([1,0,1,-11276810818409,14959107699948354572]);
assert jInvariant(E1) eq -257^3*811^3/(2^2*3^12*5^4*11);
assert jInvariant(E2) eq -441851363^3/(2^5*3*5*11*23^13);
assert jj eq {jInvariant(E1),jInvariant(E2)};
assert Conductor(E1) eq 2*3*5*11*107^2;
assert Conductor(E2) eq 23*Conductor(E1);
score := #[x : x in Eltseq(aa(E1) - aa(E2)) | x ne 0];
assert score le 1;

print "Example 8.4";

r := 1/3;
s := -3/2;
g2,h2 := Z13Equations(2,r,s);
assert IsSquare(h^2 + 4*g);
A,B,D := Z13ModuliInterpretation(2,r,s);
jprod := A^3/D;
jsum := (jprod + 1728^2 - B^2/D)/1728;
jj := {r[1] : r in Roots(X^2 - jsum*X + jprod)};
E1 := EllipticCurve([1,0,0,3796,-685680]);
E2 := EllipticCurve([1,0,0,8246713256941,11003401358367836019]);
assert jInvariant(E1) eq 29^3*61^3*103/(2^19*3^7*17);
assert jInvariant(E2) eq 13^3*103^2*270539^3/(2*3^3*17^2*19^13);
assert jj eq {jInvariant(E1),jInvariant(E2)};
assert Conductor(E1) eq 2*3*17*103^2;
assert Conductor(E2) eq 19*Conductor(E1);
score := #[x : x in Eltseq(aa(E1) - aa(E2)) | x ne 0];
assert score le 1;

print "\nThe curves of genus 0 and 1 listed in Section 2.2\n";

R<r,s> := PolynomialRing(Rationals(),2);
P<t> := PolynomialRing(Rationals());

curves := [
[
   <"P1-P1", "j = infty", s, [ t, 0, 1 ]>,
   <"P1-P1", "j = infty", r + s - 1, [ t, 1, t + 1 ]>,
   <"P1", "j = infty", r^2 + s - 1, [ t, t^2 - 1, -1 ]>,
   <"P1", "X_0(10)", r, [ 0, t, 1 ]>,
   <"P1", "X_0(25)", r + s, [ t, -t, 1 ]>,
   <"P1-P1", "", r^5 + r^4*s - 3*r^4 - r^3*s + 2*r^2*s^2 - 4*r^2*s -
        2*r*s^2 + s^3 - 4*s^2, [t*(t^2 - t + 1)*(t^2 + 2),
	-(t - 1)*(t^2 + 2)^2, (t^2 - t + 1)^2 ]>,
   <"27a1", "X_0(27)", r^2 + s, [ t, t^2, -1 ]>,
   <"36a1", "X_0(36)", r^2 + r*s - r - s + 1,
        [ t + 1, -t^2 - t - 1, -t^2 - t ]>,
   <"49a1", "X_0(49)", r^2 + r*s - s, [ t, -t^2, t - 1 ]>
],
[
   <"P1-P1", "j = infty", r, [ 0, t, 1 ]>,
   <"P1-P1", "j = infty", r^2 + r*s + r + 1, [ t^2, -t^2 - t - 1, t ]>,
   <"P1", "X_0(18)", s, [ t, 0, 1 ]>,
   <"P1-P1", "", r^2*s + r*s^2 + r*s + 2*s^2 - 2*s + 1,
        [ t^3 - t^2 + t - 1, 1, -t^2 + t ]>,
   <"19a1", "X_0(19)", r + 1, [ -1, t, 1 ]>,
   <"20a1", "X_0(20)", s - 1, [ t, 1, 1 ]>,
   <"21a1", "X_0(21)", r*s + 1, [ 1, -t^2, t ]>,
   <"24a1", "X_0(24)", r^2*s + 2*r*s + 1, [ t^3 + 2*t^2, -1, t^2 + 2*t ]>,
   <"32a1", "X_0(32)", r^2*s + r*s + 1, [ t + 1, -t^3, t^2 + t ]>,
   <"267632f1", "", r^3*s + r^2*s^2 + 3*r^2*s + 4*r*s + r + 2,
        [(t^2 + 2*t - 1)^2, -t^3*(t + 1), -t*(t^2 + 2*t - 1)]>
]];

R<x,y,z> := Universe(F);
D := [D1,D2];

function MySquareFree(poly)
  ff,c := Factorization(poly);
  assert poly eq c*(&*[P|f[1]^f[2] : f in ff]);
  c := PowerFreePart(Rationals()!c,2);
  return c*(&*[P|f[1] : f in ff | IsOdd(f[2])]);
end function;

for k in [1,2] do
  for cc in curves[k] do
    ref,lab,f0,param := Explode(cc);
    f := R!(z^Degree(f0)*Evaluate(f0,[x/z,y/z]));
    printf "%o ->",f0;
    assert Evaluate(f,param) eq 0;
    h := MySquareFree(Evaluate(F[k],param));
    if ref eq "P1-P1" then
      assert h eq 1;
    elif ref eq "P1" then
      assert Degree(h) eq 2 and LeadingCoefficient(h) eq 1;
    else
      assert Degree(h) eq 4;
      C := GenusOneModel(HyperellipticCurve(h));
      assert ref eq CremonaReference(Jacobian(C));
    end if;
    if exists(dd){ dd : dd in Table1 | dd[3] eq param
        and LegendreSymbol(k*dd[1],13) eq 1 } then
      assert lab eq Sprintf("X_0(%o)",dd[1]);
    elif GCD(D[k],f) ne 1 then
      assert lab eq "j = infty";
    else
      assert lab eq "";
    end if;
    printf " %o %o\n",ref,lab;
  end for;
end for;

print "\nTables 4 and 5";

Table4 := [[ 2, 1, -3 ], [ 3, -4, -3 ], [ 2, -5, -4 ], [ 4, 5, 3 ],
    [ 3, 1, -6 ], [ 6, -11, -9 ], [ 12, -16, -9 ], [ 14, 4, -21 ], [
    5, -23, -6 ], [ 15, -25, -18 ], [ 38, 7, 12 ], [ 15, -40, 9 ], [
    15, -63, 4 ], [ 60, -65, 16 ], [ 30, 68, 63 ], [ 21, -1, 98 ], [
    51, -136, -111 ], [ 78, -172, -169 ], [ 39, 169, 180 ], [ 56,
    -256, -245 ], [ 17, -289, -20 ], [ 95, -7, 418 ], [ 455, 169, 294
    ], [ 476, 289, 240 ], [ 104, -481, 64 ], [ 680, -175, 289 ], [
    630, 685, -324 ], [ 440, -48, 1085 ], [ 495, 81, -1144 ], [ 442,
    1224, 1183 ], [ 1430, -1469, 225 ], [ 280, -1656, 49 ], [ 1326,
    -2312, -1521 ], [ 3540, -3481, 144 ], [ 1309, -3757, 588 ], [
    4144, -999, 2695 ], [ 5110, -5329, -1176 ], [ 6552, -1352, 2835 ],
    [ 6920, -8477, -4800 ], [ 2470, 9025, 7436 ], [ 4389, -9386, -9702
    ], [ 7259, 1525, 9996 ], [ 3105, 13225, 994 ], [ 13340, 4205,
    -5819 ], [ 10540, 289, -26908 ], [ 34086, -34385, -3249 ], [ 8015,
    58166, -833 ], [ 220836, -913936, -859705 ]];

assert forall{P : P in Table4 | IsSquare(Evaluate(F[1],P))};

Table5 := [[ 1, 6, 2 ], [ 8, -5, -4 ], [ 1, -8, 6 ], [ 8, -1, 4 ], [
    2, -9, 6 ], [ 1, -7, 9 ], [ 11, -8, 22 ], [ 5, -24, 14 ], [ 27,
    -1, 12 ], [ 8, -27, -12 ], [ 20, -27, -30 ], [ 9, -40, -30 ], [
    40, -27, -24 ], [ 17, -56, 34 ], [ 15, -64, 20 ], [ 49, -64, -140
    ], [ 175, 32, 140 ], [ 64, -343, -280 ], [ 153, 343, -105 ], [
    363, -250, -165 ], [ 790, -343, 1106 ], [ 1107, -824, -246 ], [
    1176, -1331, -231 ], [ 1445, -216, 510 ], [ 532, -2197, -1235 ], [
    63, 2560, 120 ], [ 1989, -2744, -2730 ], [ 2975, 1, -255 ], [ 925,
    -2662, 4070 ], [ 175, -4608, 3080 ], [ 1007, -4913, 2584 ], [ 845,
    -4968, -1482 ], [ 5635, -6859, -1995 ], [ 7546, 1350, -735 ], [
    1682, -1331, -11484 ], [ 4563, 12167, 13455 ], [ 14175, -1331,
    4389 ], [ 1156, -15625, 5525 ], [ 13248, -42875, 21840 ], [ 78352,
    -54925, -21580 ], [ 25205, -98304, 47712 ], [ 159367, -109744,
    -81016 ]];

assert forall{P : P in Table5 | IsSquare(Evaluate(F[2],P))};

// Checking the points in Tables 4 and 5 do not lie on the
// curves of genus 0 and 1 in Section 2.2.

curves1 := [z] cat [R!(z^Degree(f[3])*Evaluate(f[3],[x/z,y/z]))
                           : f in curves[1]];
curves2 := [z] cat [R!(z^Degree(f[3])*Evaluate(f[3],[x/z,y/z]))
                           : f in curves[2]];

for P in Table4 do
  assert not exists{f : f in curves1 | Evaluate(f,P) eq 0};
end for;

for P in Table5 do
  assert not exists{f : f in curves2 | Evaluate(f,P) eq 0};
end for;

// Example 8.5

K<t> := FunctionField(Rationals());
p := t^2 - 3*t - 1;
q := t^3 - 2*t^2 - 3*t - 8;
f1 := t^7 + 4*t^6 + 8*t^5 + 6*t^4 - 8*t^3 - 24*t^2 - 27*t - 8;
f2 := t^11 + 4*t^10 + 5*t^9 - 12*t^8 - 63*t^7 - 124*t^6
    - 137*t^5 - 80*t^4 - 61*t^3 - 72*t^2 - 153*t + 8;
g1 := t^31 + 23*t^30 + 270*t^29 + 2379*t^28 + 17607*t^27 + 110676*t^26
  + 586710*t^25 + 2624262*t^24 + 9977316*t^23 + 32555542*t^22 +
  92002244*t^21 + 226872066*t^20 + 490871649*t^19 + 935166681*t^18 +
  1571157252*t^17 + 2326844467*t^16 + 3029704865*t^15 +
  3450459162*t^14 + 3407984048*t^13 + 2880044002*t^12 +
  2037108963*t^11 + 1159162859*t^10 + 486247810*t^9 + 109783239*t^8 -
  25731445*t^7 - 37205624*t^6 - 17036352*t^5 - 3782272*t^4 - 99968*t^3
  + 90624*t^2 + 50176*t + 8192;
g2 := t^46 + 34*t^45 + 586*t^44 + 6298*t^43 + 44058*t^42 +
  176260*t^41 - 109414*t^40 - 8297186*t^39 - 80245684*t^38 -
  525645970*t^37 - 2745467763*t^36 - 12123387334*t^35 -
  46607587414*t^34 - 158717868270*t^33 - 484152170146*t^32 -
  1332985962708*t^31 - 3330160262195*t^30 - 7577440618404*t^29 -
  15743853949526*t^28 - 29918470706666*t^27 - 52043929525145*t^26 -
  82879264209560*t^25 - 120751034770587*t^24 - 160726419562280*t^23 -
  194995600204665*t^22 - 214887899363332*t^21 - 214056270394992*t^20 -
  191412686703724*t^19 - 152125444758091*t^18 - 105839083071274*t^17 -
  62871218212945*t^16 - 30397468667742*t^15 - 10590883463569*t^14 -
  1344855906720*t^13 + 1389594877516*t^12 + 1310083979082*t^11 +
  634627697211*t^10 + 183457593768*t^9 + 16079661600*t^8 -
  12752378496*t^7 - 6959066304*t^6 - 1639051008*t^5 - 144015872*t^4 -
  11264000*t^3 - 10680320*t^2 - 2752512*t - 262144;
E1 := EllipticCurve([-3*p^2*q*f1,2*p^2*q^2*f2]);
E2 := EllipticCurve([-3*p^3*q^3*g1,2*p^4*q^5*g2]);
jjdir := [jInvariant(E1),jInvariant(E2)];

// Example 8.6

p := t^2 + t + 1;
q := 5*t^2 + 8*t + 11;
f1 := t^4 - 13*t^3 - 4*t^2 - 5*t + 1; 
f2 := 59*t^9 + 183*t^8 + 477*t^7 + 315*t^6 + 54*t^5
     - 570*t^4 - 499*t^3 - 429*t^2 - 123*t - 43;
g1 := 211*t^14 + 665*t^13 + 1079*t^12 + 414*t^11 - 1754*t^10 -
    5658*t^9 - 9756*t^8 - 12536*t^7 - 12796*t^6 - 10606*t^5
    - 7358*t^4 - 4030*t^3 - 1831*t^2 - 553*t - 131;
g2 := 3107*t^23 + 45563*t^22 + 257591*t^21 + 945191*t^20 +
    2573711*t^19 + 5625383*t^18 + 10317391*t^17 + 16362175*t^16 +
    22933723*t^15 + 28854995*t^14 + 32992127*t^13 + 34527311*t^12 +
    33211633*t^11 + 29341897*t^10 + 23725261*t^9 + 17433613*t^8 +
    11528413*t^7 + 6782485*t^6 + 3490781*t^5 + 1545437*t^4 +
    568361*t^3 + 168353*t^2 + 35789*t + 4973;
E1 := EllipticCurve([-3*p*q^2*f1,2*p*q^2*f2]);
E2 := EllipticCurve([3*p*q^3*g1,2*p*q^4*g2]);
jjskew := [jInvariant(E1),jInvariant(E2)];

print "\nTable 6 : Pairs of 13-congruent elliptic curves\n";

// Each row lists: k, ref1, ref2, W1, W2, P, label

// k = 1 or 2 according as the congruence is direct or skew
// ref1, ref2 : Cremona labels for E1 and E2
// W1, W2 : coefficients of a Weierstrass equation for E1 and E2
// P : the corresponding point on Z(13,k) 
// label : some information about how the examples was found
//         (either a value of t, for examples coming from the
//          P^1-families in Examples 8.5 and 8.6, or "Table 4/5"
//	   if P comes from one of those tables, or "ec - rank 2"
//	   for examples coming from the family parametrised by
//	   a rank 2 elliptic curve (see Section 2.2))

Table6 := [
<2, "52a1", "988b1", [ 0, 0, 0, 1, -10 ], [ 0, 0, 0, -362249, 165197113 ], [ 8, -1, 4 ], "Table 5">,
<1, "52a2", "988b1", [ 0, 0, 0, -4, -3 ], [ 0, 0, 0, -362249, 165197113 ], [ 15, -25, -18 ], "Table 4">,
<1, "345b1", "10005m1", [ 0, 1, 1, -1, 1 ], [ 0, 1, 1, -104207741, 135268278965 ], [ 12, -16, -9 ], "Table 4">,
<2, "735c1", "9555h1", [ 0, -1, 1, 5, 6 ], [ 0, -1, 1, 654855, -11895892 ], [ 1, -8, 6 ], "Table 5">,
<2, "735c2", "9555h1", [ 0, -1, 1, -205, 1203 ], [ 0, -1, 1, 654855, -11895892 ], [ 15, -64, 20 ], "Table 5">,
<1, "1190a1", "265370d1", [ 1, -1, 0, -145, -2979 ], [ 1, -1, 0, -4246061054705, 4353350273721622301 ], [ 39, 169, 180 ], "Table 4">,
<2, "1190a2", "265370d1", [ 1, -1, 0, -4145, -101379 ], [ 1, -1, 0, -4246061054705, 4353350273721622301 ], [ 175, 32, 140 ], "Table 5">,
<1, "1274h1", "21658t1", [ 1, 0, 0, -2024, 34880 ], [ 1, 0, 0, -118742534, -5724570079852 ], [ 21, -1, 98 ], "Table 4">,
<1, "1445b1", "10115e1", [ 1, 1, 0, 717, 5422 ], [ 1, 1, 0, -506478, 86623057 ], [ 2, -5, -4 ], "Table 4">,
<2, "1960i1", "21560l1", [ 0, 1, 0, -11776, 568224 ], [ 0, 1, 0, 5764344, -2909997035 ], [ 20, -27, -6 ], "t = -3">,
<1, "3990ba1", "43890cu1", [ 1, 0, 0, -190366575, 325694589866937 ], [ 1, 0, 0, 5080055, 9229629437 ], [ 476, 289, 240 ], "Table 4">,
<2, "4719b1", "33033k1", [ 0, -1, 1, -394137, -95108641 ], [ 0, -1, 1, -72034607, 235356539186 ], [ 8, -27, -12 ], "Table 5">,
<2, "5070j1", "35490bg1", [ 1, 0, 1, -836554, -294572068 ], [ 1, 0, 1, -31671635129, 2169468305732756 ], [ 27, -1, 12 ], "Table 5">,
<1, "11638o1", "151294h1", [ 1, -1, 1, -1951, 41799 ], [ 1, -1, 1, -17498345946, -890914836928711 ], [ 6, -11, -9 ], "Table 4">,
<2, "12274c1", "135014s1", [ 1, 1, 0, -26, 46 ], [ 1, 1, 0, -11932483101, -493031790814211 ], [ 2, -686, -21 ], "ec - rank 2">,
<1, "13688b1", "8363368*", [ 0, -1, 0, -1152, 15212 ], [ 0, -1, 0, -25715622098537, -29022548679301803715 ], [ 315, -324, 49 ], "t = 4">,
<2, "14175k1", "184275o1", [ 1, -1, 1, -14555, 691822 ], [ 1, -1, 1, 28536070, -44642935678 ], [ 5, -8, -2 ], "t = -2">,
<1, "20184i1", "20184j1", [ 0, 1, 0, -162593, -25306749 ], [ 0, 1, 0, -26212568, 51646275696 ], [ 3, -4, -3 ], "Table 4">,
<2, "23660f1", "733460*", [ 0, -1, 0, 5859, 792466 ], [ 0, -1, 0, -90590037581, -10631522405253959 ], [ 40, -1, 12 ], "t = 1/3">,
<1, "27930s1", "27930r1", [ 1, 1, 0, 38, -146 ], [ 1, 1, 0, -3695905122277, 2735775223839669901 ], [ 231, -242, 49 ], "t = 2">,
<2, "29970f1", "4705290*", [ 1, -1, 0, -5724, -166832 ], [ 1, -1, 0, -3077249264559, 2071858234462281613 ], [ 64, -343, -280 ], "Table 5">,
<2, "69230m1", "761530*", [ 1, 1, 0, -14562, 106188916 ], [ 1, 1, 0, -237573527324147, -8269595392408362433691 ], [ 5635, -6859, -1995 ], "Table 5">,
<1, "80408l1", "8282024*", [ 0, -1, 0, -12320, -521876 ], [ 0, -1, 0, 25744150157415, -21222569268437349971 ], [ 429, -484, -169 ], "t = -4">,
<1, "83030b1", "913330*", [ 1, 1, 0, 12, -6832 ], [ 1, 1, 0, -93645356988, 11030012229735752 ], [ 3, 9, 2 ], "t = -1/2">,
<2, "92950q1", "2881450*", [ 1, 0, 1, -2431576, 1459354548 ], [ 1, 0, 1, -6594535339396, -6517800906146351742 ], [ 25, 27, -12 ], "t = -3/4">,
<1, "95370cl1", "95370cm1", [ 1, 1, 1, -7941596025, 272489201594055 ], [ 1, 1, 1, -9560, -363763 ], [ 3, 1, -6 ], "Table 4">,
<2, "98010s1", "98010t1", [ 1, -1, 0, -4966044, -4304874700 ], [ 1, -1, 0, -2288259, 1332884053 ], [ 5, 1, -2 ], "t = -1/2">,
<2, "162266e1", "17686994*", [ 1, 0, 1, -91, -186 ], [ 1, 0, 1, 2014602555250579, 13649307693566295203064 ], [ 1, 6, 2 ], "Table 5">,
<2, "177735a1", "1955085*", [ 1, 0, 0, -7596973640, -254864925743775 ], [ 1, 0, 0, -8646107309665, 9785409709784043350 ], [ 153, 343, -105 ], "Table 5">,
<2, "185900a1", "7621900*", [ 0, 1, 0, 91542, 23170213 ], [ 0, 1, 0, 3428811305292, 9985577786046319588 ], [ 40, -27, 12 ], "t = 3">,
<2, "237538k1", "23041186*", [ 1, 0, 1, -144655596, -669667246038 ], [ 1, 0, 1, -1144586636657031, 14904885221715681198846 ], [ 2209, -6655, -3102 ], "ec - rank 2">,
<2, "255162e1", "4848078*", [ 1, 1, 0, -110978, -33359766 ], [ 1, 1, 0, 1196689182648437, 487962218257288589629 ], [ 300, -343, -42 ], "t = -7">,
<1, "271545f1", "589524195*", [ 0, 1, 1, -60, 151229 ], [ 0, 1, 1, -85225672741555130, -9750594868956521658113341 ], [ 308, 484, 147 ], "t = -1/3">,
<2, "314330i1", "314330j1", [ 1, 1, 0, -2060727, 1137761789 ], [ 1, 1, 0, -182683503892, -30053755938589104 ], [ 40, -27, -24 ], "Table 5">,
<2, "1082118*", "20560242*", [ 1, 0, 0, 3796, -685680 ], [ 1, 0, 0, 8246713256941, 11003401358367836019 ], [ 2, -9, 6 ], "Table 5">,
<2, "1137150*", "76189050*", [ 1, -1, 1, 51175, -27858823 ], [ 1, -1, 1, -18471642425, 1412242059112927 ], [ 5, -24, 14 ], "Table 5">,
<2, "1242150*", "1242150*", [ 1, 1, 1, 1602, 42531 ], [ 1, 1, 1, -2305248, -1433221119 ], [ 15, -1, 6 ], "t = 1/2">,
<1, "1296924*", "1437132*", [ 0, 1, 0, -19888964357, 1068009894095943 ], [ 0, 1, 0, 9417106468, -634288021724460 ], [ 5, -23, -6 ], "Table 4">,
<2, "1425720*", "1425720*", [ 0, 1, 0, -307471160, 2075072352048 ], [ 0, 1, 0, -1601087520, -25155895637007 ], [ 20, -27, -30 ], "Table 5">,
<2, "1695978*", "429082434*", [ 1, -1, 1, -11, -1061 ], [ 1, -1, 1, -178678951443101, 920087202239659475983 ], [ 400, -1, 56 ], "t = 1/7">,
<2, "2141594*", "49256662*", [ 1, 1, 1, -189988, 31824757 ], [ 1, 1, 1, 4050708485382, -1612141371345940345 ], [ 13, -27, -6 ], "t = -3/2">,
<2, "2147950*", "2147950*", [ 1, 0, 1, -8501, 301638 ], [ 1, 0, 1, -945422025951, -353823578680306702 ], [ 68, 27, -30 ], "t = -3/5">,
<2, "2164218*", "413365638*", [ 1, 1, 0, -18552, -992448 ], [ 1, 1, 0, -11741563342271397922, 15774354781507281862095272422 ], [ 1, -7, 9 ], "Table 5">,
<2, "2228037*", "69069147*", [ 0, -1, 1, -189547, -32282937 ], [ 0, -1, 1, -64962009727347, 202002358618269858479 ], [ 17, -56, 34 ], "Table 5">,
<1, "2428110*", "31565430*", [ 1, -1, 0, -12807848685, -557903432180059 ], [ 1, -1, 0, 1971001530, 58324261594320 ], [ 56, -256, -245 ], "Table 4">,
<2, "3647770*", "69307630*", [ 1, 0, 0, 188184529130, -15585854275904988 ], [ 1, 0, 0, -620002347510, -252303438377376890 ], [175,-4608,3080],"Table 5">,
<1, "3778170*", "86897910*", [ 1, 0, 1, -464619, -122105558 ], [ 1, 0, 1, -11276810818409, 14959107699948354572 ], [ 4, 5, 3 ], "Table 4">,
<1, "3944850*", "3944850*", [ 1, 1, 0, -64501465, -212536994075 ], [ 1, 1, 0, -23675323765, 3548135402525425 ], [ 60, -65, 16 ], "Table 4">,
<2, "4083510*", "730948290*", [ 1, 1, 0, -296544158, -1965664860858 ], [ 1, 1, 0, -195637371732362563, 33311913622399107383507917 ], [ 790, -343, 1106 ], "Table 5">,
<2, "4746924*", "507920868*", [ 0, 0, 0, -8568, -523404 ], [ 0, 0, 0, 1048664270367, 82155952824381009 ], [ 104, -125, -20 ], "t = -5">,
<1, "5429670*", "320350530*", [ 1, 1, 0, -5857, -176459 ], [ 1, 1, 0, -105056013817601457, 13105647232667475016370301 ], [ 357, -289, 98 ], "t = 1/2">,
<2, "7495800*", "397277400*", [ 0, 1, 0, 10592, -19622512 ], [ 0, 1, 0, -699167391833, 225520483266015963 ], [ 156, -1, 30 ], "t = 1/5">,
<1, "7707798*", "27925352154*", [ 1, -1, 0, -39048, 2878784 ], [ 1, -1, 0, 19869245527052530992, -277948797400052128998533765440 ], [ 665, -361, 196 ], "t = 1/4">,
<2,"10052196*","3608738364*", [ 0, -1, 0, -10189516, 12522651352 ], [ 0, -1, 0, -8976490985592664259381, -329204023344019125753175545966726 ], [ 13248, -42875, 21840 ],"Table 5">,
<1, "15211515*", "1566786045*", [ 0, 1, 1, -43446, -3449539 ], [ 0, 1, 1, 1113776489639155544, 564644163389589788997195325 ], [ 84, -108, -49 ], "t = -3">,
<2, "15893150*", "842336950*", [ 1, 0, 0, -2503259978878, 1522913403637913412 ], [ 1, 0, 0, 611451736237, -171370262183016233 ], [ 4563, 12167, 13455 ],"Table 5">,
<2, "16207776*", "1183167648*", [ 0, 0, 0, -528723, -148009914 ], [ 0, 0, 0, -11693822131368, 15391376935061817456 ], [ 148, 125, -70 ], "t = -5/7">,
<2, "17859765*", "553652715*", [ 1, 0, 1, 36402, 2596003 ], [ 1, 0, 1, 61474275312, 345919051351921891 ], [ 15, -8, 6 ], "t = 2">,
<1, "21140427*", "264234197073*", [ 1, 0, 0, -13696, -625063 ], [ 1, 0, 0, -3618683202024979848780461, -5993415420929929580089303632892835256 ], [ 38, 7, 12 ], "Table 4">,
<1, "21219400*", "233413400*", [ 0, 0, 0, -1735235, -879803250 ], [ 0, 0, 0, -2128342849319635, 37792924308868253391975 ], [ 15, -40, 9 ], "Table 4">,
<2, "21997290*", "285964770*", [ 1, 1, 0, -11232482, 226369891914 ], [ 1, 1, 0, 3090316166263883, -100105541160723716277299 ], [1176,-1331,-231],"Table 5">,
<2, "22817366*", "22817366*", [ 1, 1, 0, -89573, -10391909 ], [ 1, 1, 0, 113451387, 311217061069 ], [ 68, -125, -30 ], "t = -5/3">,
<2, "27674878*", "5451950966*", [ 1, 0, 1, -5492, -156686 ], [ 1, 0, 1, -38717478800211413318757, 2932154316553796384854275640602366 ], [ 2975, 1, -255 ], "Table 5">,
<2, "28958240*", "8658513760*", [ 0, -1, 0, 358239, -110361935 ], [ 0, -1, 0, 4101633879379, -36050439858440458139 ], [ 580, -343, 210 ], "t = 7/3">,
<1, "29402282*", "2087562022*", [ 1, 0, 1, -1304255, 573204114 ], [ 1, 0, 1, -13749415080437235, -659405306400652793219314 ], [ 15, -63, 4 ], "Table 4">,
<2, "31744713*", "2507832327*", [ 1, 0, 0, 21959, -7753228 ], [ 1, 0, 0, 1972761107584, -1765374812094881943 ], [ 156, -125, 30 ], "t = 5">,
<2, "32301984*", "8107797984*", [ 0, -1, 0, -21504968, 199794217908 ], [ 0, -1, 0, 12528698835589907, 140361352934828532115621 ], [ 6516, -6859, -342 ], "t = -19">,
<2, "35322287*", "671123453*", [ 1, 0, 0, -393303, 94966976 ], [ 1, 0, 0, -22852599268038, -42048542271993472027 ], [ 13, 8, -6 ], "t = -2/3">,
<1, "43998465*", "43998465*", [ 0, 1, 1, -11815876101, 494360253359921 ], [ 0, 1, 1, -1045160614115401, -13005382605683716784954 ], [ 30, 68, 63 ], "Table 4">,
<2, "45269037*", "316883259*", [ 1, -1, 0, -343529853, -2450640574576 ], [ 1, -1, 0, -125442562938, 17210082757192829 ], [ 11, -8, 22 ], "Table 5">,
<2, "45349590*", "1315138110*", [ 1, 0, 1, -4263559, -4105166518 ], [ 1, 0, 1, -38288679841537139, 2871254923133319032617712 ], [ 9, -40, -30 ], "Table 5">,
<1, "46427580*", "5448415795740*", [ 0, 0, 0, 87588, 112263309 ], [ 0, 0, 0, -155167893413948529579310497, -743968050812936330060724078112954681011 ], [ 40807, -42632, 9245 ], "t = 8/5">,
<1, "48963840*", "42941287680*", [ 0, 1, 0, 1169, -11635 ], [ 0, 1, 0, -4499864313926385501, -4041445170634540074426378501 ], [ 936, -972, 169 ], "t = 3">,
<2, "53420780*", "101018694980*", [ 0, 1, 0, -31565, -2170225 ], [ 0, 1, 0, 118726245645086070, 6881780210595535740908405 ], [ 1000, -2197, -468 ], "t = -13/9">,
<2, "56756700*", "1645944300*", [ 0, 0, 0, -29925, -2710575 ], [ 0, 0, 0, -158749452642375, 769868361733048278750 ], [ 104, 1, -20 ], "t = -1/5">,
<2, "61455779*", "22062624661*", [ 1, 1, 0, 12219, -132986 ], [ 1, 1, 0, -1804327287789336, -151148099350530761089211 ], [ 175, -64, 84 ], "t = 4/3">,
<1, "62310225*", "2554719225*", [ 0, 1, 1, -235359433, -1389856759181 ], [ 0, 1, 1, -524281576902799133, 156023576520368819246642594 ], [ 78, -172, -169 ], "Table 4"> ];

function IsInP1Family(E1,E2,k)
  jj := k eq 1 select jjdir else jjskew;
  for m in [1,2] do
    rts := Roots(Numerator(jj[1]-jInvariant(E1)));
    if #rts ne 0 then
      assert #rts eq 1;
      t := rts[1,1];
      assert Evaluate(jj[1],t) eq jInvariant(E1);
      assert Evaluate(jj[2],t) eq jInvariant(E2);
      return true,t;
    end if;
    E1,E2 := Explode([E2,E1]);
  end for;
  return false,_;
end function;

R<x,y,z> := PolynomialRing(Rationals(),3);
P2 := Proj(R);
g1,h1 := Z13Equations(1,x/z,y/z);
F1 := R!(z^12*(4*g1 + h1^2));
g2,h2 := Z13Equations(2,x/z,y/z);
F2 := R!(z^14*(4*g2 + h2^2));
F := [F1,F2];

A1,B1,D1 := Z13ModuliInterpretation(1,x/z,y/z);
jprod1 := A1^3/D1;
jsum1 := (jprod1 + 1728^2 - B1^2/D1)/1728;
A2,B2,D2 := Z13ModuliInterpretation(2,x/z,y/z);
jprod2 := A2^3/D2;
jsum2 := (jprod2 + 1728^2 - B2^2/D2)/1728;
jsum := [jsum1,jsum2];
jprod := [jprod1,jprod2];

// The equations for the curves in P^2 whose inverse images
// in Z(13,k) give the known infinite families of pairs of
// 13-congruent elliptic curves

polys := [ x^5 + x^4*y - 3*x^4*z - x^3*y*z + 2*x^2*y^2*z - 4*x^2*y*z^2
    - 2*x*y^2*z^2 + y^3*z^2 - 4*y^2*z^3, x^2*y + x*y^2 + x*y*z +
    2*y^2*z - 2*y*z^2 + z^3, x^3*y + x^2*y^2 + 3*x^2*y*z + 4*x*y*z^2 +
    x*z^3 + 2*z^4 ];

P<X> := PolynomialRing(Rationals());
jseen := [];
toprint := [];

for ctr in [1..#Table6] do

  dd := Table6[ctr];
  k,ref1,ref2,W1,W2,pt,label := Explode(dd);
  E1 := EllipticCurve(W1);
  E2 := EllipticCurve(W2);
  assert CR(E1) eq ref1;
  assert CR(E2) eq ref2;
  printf "%o %o ",CR(E1),CR(E2);

// Checking some traces of Frobenius

  diffr := aa(E1) - aa(E2);
  assert #[x : x in Eltseq(diffr) | x ne 0] lt 5;

// Checking the example is new

  assert jInvariant(E1) notin jseen;
  if jInvariant(E2) in jseen then
    assert (ctr eq 2 and CR(E2) eq "988b1")
        or (ctr eq 5 and CR(E2) eq "9555h1")
        or (ctr eq 7 and CR(E2) eq "265370d1");
  end if;
  assert #IsogenousCurves(E1) eq 1 xor ctr in [1,2,4,5,6,7];
  assert #IsogenousCurves(E2) eq 1;
  assert #TorsionSubgroup(E1) eq 1 xor ctr in [1,2,6,7];
  assert #TorsionSubgroup(E2) eq 1;
  assert Conductor(MinimalTwist(E1)) eq Conductor(E1);
  assert Conductor(MinimalTwist(E2)) ge Conductor(E1);
  assert not IsIsogenous(E1,E2);
  dp := DivisionPolynomial(E1,13);
  assert IsIrreducible(dp);
  assert jInvariant(E1) ne jInvariant(E2);
  Append(~jseen,jInvariant(E1));
  Append(~jseen,jInvariant(E2));

// Checking E1 and E2 correspond to the given point on Z(13,k)

  assert IsSquare(Evaluate(F[k],pt));
  jpoly := X^2 - Evaluate(jsum[k],pt)*X + Evaluate(jprod[k],pt);
  rts := Roots(jpoly);
  assert {r[1] : r in rts} eq {jInvariant(E1),jInvariant(E2)};
  assert GCD(pt) eq 1;

// Checking the description of how the example was found

  flag,t := IsInP1Family(E1,E2,k);
  rr := [i : i in [1..3] | Evaluate(polys[i],pt) eq 0];

  if flag then
    assert label eq Sprintf("t = %o",t);
    assert rr eq [k];
  elif label eq "Table 4" then
    assert k eq 1 and pt in Table4 and rr eq [];
  elif label eq "Table 5" then
    assert k eq 2 and pt in Table5 and rr eq [];
  else
    assert label eq "ec - rank 2";
    assert rr eq [3];
  end if;
  print label;

// Printing out the .tex file version of Table 6

  if CR(E1) notin {"52a2","735c2","1190a2"} then 
    if not flag then 
      Append(~toprint,Sprintf("%o & %o & %o ",k,CR(E1),CR(E2)));
    else
      Append(~toprint,Sprintf("%o & {\\bf{%o}} & {\\bf{%o}} ",k,CR(E1),CR(E2)));
    end if;
  end if;

//  if ctr mod 20 eq 0 then print "continue?"; read dit; end if;

end for;

print "";

ncols := 29;
for m in [1..ncols] do
  printf "%o & %o \\\\ \n",toprint[m],toprint[ncols + m];
end for;

print "\nAppendix A\n";

K<t> := FunctionField(Rationals());

// The j-map on X_0(13) (following Fricke) 
j13 := (t^2 + 5*t + 13)*(t^4 + 7*t^3 + 20*t^2 + 19*t + 1)^3/t;

// The genus 2 curve X_1(13) and its automorphisms
P<x> := PolynomialRing(Rationals());
h := (x^3 - 3*x + 1)^2 - 2*(x^3 - 3*x + 1)*x*(x - 1) + 5*x^2*(x - 1)^2;
assert Evaluate(h,1/(1-x)) eq h/(1-x)^6;
X1 := HyperellipticCurve(h);
assert #Automorphisms(X1) eq 6;

// The surface Sigma in the proof of Theorem A.1
R<T,X> := PolynomialRing(Rationals(),2);
la := T^3 - 3*T + 1;
mu := T*(T - 1);
f := (X^2 - 2*X + 5)*((la^2 - 2*la*mu + 5*mu^2)*X^2
  - 2*(la^2 + la*mu + 6*mu^2)*X + (5*la^2 - 12*la*mu + 72*mu^2));

K1<x1,y1> := FunctionField(X1);
X1K1 := BaseChange(X1,K1);
K2<x2,y2> := FunctionField(X1K1);

// Formulae for the map alpha
T := (x1*x2 - x1 + 1)/(x2 - x1);
X := (x2^3 - 3*x2 + 1)/(x2*(x2 - 1));
Y := (X^2 - 3*X + 9)*y1*y2/(x1 - x2)^3;
assert Y^2 eq Evaluate(f,[T,X]);

// Formulae for the map beta
s1 := ((T^3 - 3*T + 1)*X - 9*T*(T - 1))
     /(T*(T - 1)*X + (T^3 - 3*T^2 + 1));
s2 := X;

// Checking the diagram commutes
assert s1 eq (x1^3 - 3*x1 + 1)/(x1*(x1 - 1));
assert s2 eq (x2^3 - 3*x2 + 1)/(x2*(x2 - 1));

// Exhibiting a fibre of positive rank
T := -2;
E := EllipticCurve([0,-1,0,-685532,-213343020]);
P1 := E![ 11523732043512033/2408363^2, 1105503446246122991427360/2408363^3 ];
KE<x,y> := FunctionField(E);
X := (1768645*y + 30366*x^2 + 4012598*x - 6595451308)
          /(15876*x^2 - 10801141*x - 7669747486);
assert IsSquare(Evaluate(f,[T,X]));

// Example A.3
T := 17/33;
X := 1;
Y := 126340/33^3;
assert Y^2 eq Evaluate(f,[T,X]);
s1 := ((T^3 - 3*T + 1)*X - 9*T*(T - 1))
     /(T*(T - 1)*X + (T^3 - 3*T^2 + 1));
s2 := X;
E1 := EllipticCurve([0,-1,1,-2,-1]);
E2 := EllipticCurve([0,-1,1,-1424883795842044404862,
                         -20702237422068075268318817670099]);
assert Evaluate(j13,s1 - 4) eq jInvariant(E2);
assert Evaluate(j13,s2 - 4) eq jInvariant(E1);
print E1;
print E2;

print "";

