
// 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 curves X_E(13,1) and X_E(13,2)

// The calculations are checked in roughly the order of the paper

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

load "congr13-formulae.m";

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

print "\nSection 2.1\n";

function ApplyDiff(f,g)
// Computes f( d/dx_1, ..., d/dx_7 ) g( x_1, ..., x_7)
  ans := 0;
  for t in Terms(f) do
    dd := [i : j in [1..Degree(t,i)], i in [1..7]];
    g0 := g;
    for d in dd do
      g0 := Derivative(g0,d);
    end for;
    ans +:= Coefficients(t)[1]*g0;
  end for;
  return ans;
end function;

function DualQuadrics(ff)
// Compute a basis for the vector space
// { g : f( d/dx_1, ..., d/dx_7 ) g = 0 for all f in ff } 
  R := Universe(ff);
  k0 := BaseRing(R);
  mons := MD(R,2);
  mat := Matrix(k0,28,14,[ApplyDiff(f,m): f in ff, m in mons]);
  km := KernelMatrix(mat);
  assert Nrows(km) eq 14;
  gg := [&+[km[i,j]*mons[j]: j in [1..#mons]]: i in [1..14]];
  assert forall{[f,g] : f in ff,g in gg | ApplyDiff(f,g) eq 0};
  return gg;
end function;

function SyzygySupport(ff)
// Compute the support of all linear syzygies 
  R :=  Universe(ff);
  assert #ff eq 14;
  M := RSpace(R,14);
  bb := [R.i*M.j : i in [1..7],j in [1..14]];
  mons3 := MD(R,3);
  mat := Matrix(#bb,#mons3,[[MC(f,m): m in mons3]
         where f is &+[b[i]*ff[i] : i in [1..14]] : b in bb]);
  km := KernelMatrix(mat);
  assert Nrows(km) eq 14;
  SYZ := [&+[km[i,j]*bb[j] : j in [1..#bb]]: i in [1..14]];
  M := Transpose(Matrix(SYZ));
  assert Vector(ff)*M eq 0;
  V := VectorSpace(BaseRing(R),14);
  U := sub<V|>;
  for i in [1..14] do
    for j in [1..7] do
      vec := V![MC(M[k,i],R.j): k in [1..14]];
      U +:= sub<V|vec>;
      if Dimension(U) eq 13 then break i; end if;
    end for;
  end for;
  km := Matrix(Basis(U));  
  return [&+[km[i,j]*ff[j]: j in [1..#ff]]: i in [1..13]];
end function;  

function GetCubics(ff)
// Compute a basis for the space of cubics whose partial derivatives
// belong to the space of quadrics spanned by ff
  R := Universe(ff);
  k0 := BaseRing(R);
  mons2 := MD(R,2);
  mat := Matrix(#ff,#mons2,[MC(f,m): m in mons2,f in ff]);
  V := VectorSpace(k0,#mons2);
  U := sub<V|[V!mat[i] : i in [1..Nrows(mat)]]>;
  UU,quomap := quo<V|U>;
  tovec := func<f|Vector(#mons2,[MC(f,m): m in mons2])>;
  tovec1 := func<f|Eltseq(quomap(tovec(f)))>;
  mons3 := MD(R,3);
  mat := Matrix(#mons3,7*Dimension(UU),
     [&cat[tovec1(Deriv(m,j)): j in [1..7]]: m in mons3]);
  km := KernelMatrix(mat);
  assert Nrows(km) eq 7;
  return [&+[km[i,j]*mons3[j]: j in [1..#mons3]]: i in [1..7]];
end function;

function HaveSameSpan(forms1,forms2,d)
// Check the two sets of forms (of degree d) span the same
// vector space
  R := Universe(forms1);
  mons := MD(R,d);
  V := VectorSpace(BaseRing(R),#mons);
  W1 := sub<V|[V![MC(f,m): m in mons]: f in forms1]>;
  W2 := sub<V|[V![MC(f,m): m in mons]: f in forms2]>;
  return W1 eq W2;
end function;

function ComputeCubics(quadsU)
// Compute the space of cubics in Theorem 2.1
  assert #quadsU eq 14;
  print "  14 quadrics on (P^6)^*";
  quadsUperp := DualQuadrics(quadsU);
  assert #quadsUperp eq 14;
  print "   --> 14 quadrics on P^6";
  quadsV := SyzygySupport(quadsUperp);
  assert #quadsV eq 13;
  print "   --> 13 quadrics on P^6";
  cubics := GetCubics(quadsV);
  assert #cubics eq 7;
  print "   --> 7 cubics on P^6";
  return cubics;
end function;

print "Theorem 2.1 : Equations for X(13)";

K<z13> := CyclotomicField(13);
rr := [0,1,4,3,12,9,10];
P6 := ProjectiveSpace(Rationals(),6);
pts := {P6(K)![1,0,0,0,0,0,0]};
pts join:= {P6(K)![z13^(r*s): r in rr]: s in [0..12]};
quadsU := MinimalBasis(Cluster(pts));
cubics := ComputeCubics(quadsU);
R<x0,x1,x2,x3,x4,x5,x6> := CoordinateRing(P6);
cubics0 := [
    -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
];
assert HaveSameSpan(cubics,cubics0,3);
print "  Agrees with answer in proof of Theorem 5.2";

print "Formulae for the 14 hyperplanes";

K<t> := FunctionField(Rationals());
// The j-map on X_0(13) (following Fricke)
j := (t^2 + 5*t + 13)*(t^4 + 7*t^3 + 20*t^2 + 19*t + 1)^3/t;
c4 := (t^2 + 5*t + 13)*(t^2 + 6*t + 13)*(t^4 + 7*t^3 + 20*t^2 + 19*t + 1);
c6 := (t^2 + 5*t + 13)*(t^2 + 6*t + 13)^2
            *(t^6 + 10*t^5 + 46*t^4 + 108*t^3 + 122*t^2 + 38*t - 1);
E := EllipticCurve([-27*c4,-54*c6]);
assert jInvariant(E) eq j;

P<s,t> := PolynomialRing(Rationals(),2);
c4 := P!Evaluate(c4,t);
c6 := P!Evaluate(c6,t);

// The polynomials f1,..,f7 and g1,..,g7 defined before Theorem 2.2.

f1 := 1;
f2 := -(t + 1);
f3 := 3*s*(t + 2)*(t^2 + 5*t + 13)*(t^2 + 6*t + 13);
f4 := -9*s*(t + 1)*(t^2 + 5*t + 13)*(t^2 + 6*t + 13);
f5 := 108*s^2*(t^2 + 5*t + 13)*(t^2 + 6*t + 13)
            *(t^3 + 5*t^2 + 10*t + 2) + 27*s^2*(t + 3)*c4;
f6 := 162*s^2*(t^2 + 5*t + 13)*(t^2 + 6*t + 13)
            *(t^3 + 6*t^2 + 14*t + 7) + 27*s^2*(t + 4)*c4;
f7 := 11664*s^3*(t + 1)*(t^2 + 5*t + 13)*(t^2 + 6*t + 13)^2
            + 54*s^2*c4*f4;

g1 := 2;
g2 := 2*(t + 1);
g3 := 3*s*(t^2 + 6*t + 13)*(t^3 + 4*t^2 + 8*t - 1);
g4 := 12*s*(t^2 + 6*t + 13)*(t^2 + 3*t + 5);
g5 := 6*s*(t^2 + 6*t + 13)*(t^3 + 8*t^2 + 20*t + 7);
g6 := 108*s^2*(t + 1)^2*(t^2 + 5*t + 13)*(t^2 + 6*t + 13)
	    - 9*s^2*(t + 3)*c4;
g7 := -216*s^2*(t - 1)*(t^2 + 5*t + 13)*(t^2 + 6*t + 13)
            - 18*s^2*(t + 2)*c4;

ff := [[f1,f2,f3,f4,f5,f6,f7],[g1,g2,g3,g4,g5,g6,g7]];

A := -27*s^2*c4;
B := -54*s^3*c6;

K<a,b> := FunctionField(Rationals(),2);
E := EllipticCurve([a,b]);
P<t> := PolynomialRing(K);
poly := P!(t*(Evaluate(j,t) - jInvariant(E)));
L<t> := quo<P|poly>;
s := b/a*Evaluate(A/B,[1,t]);
assert Evaluate(A,[s,t]) eq a;
assert Evaluate(B,[s,t]) eq b;

R<[x]> := PolynomialRing(K,7);
mons := MD(R,2);

for k in [1,2] do
  printf "Theorem 2.2 : Equations for X_E(13,%o)\n",k;
  pt := [L|Evaluate(f,[s,t]): f in ff[k]];
  mat := Matrix(#mons,14,[EltseqPad(Evaluate(m,pt)) : m in mons]);
  km := KernelMatrix(mat);
  assert Nrows(km) eq 14;
  quadsU := [&+[km[i,j]*mons[j] : j in [1..#mons]]: i in [1..14]];
  cubics := ComputeCubics(quadsU);
  cubics1 := XE13Cubics(k,a,b,x);
  assert HaveSameSpan(cubics,cubics1,3);
  print "  Agrees with answer in congr13-formulae.m";
end for;

print "Remark 2.3 : Checking the gradings";
K<a,b,t> := FunctionField(Rationals(),3);
R<[x]> := PolynomialRing(K,7);
ww := [[3,3,2,2,1,1,0],[2,2,1,1,1,0,0]];
dd := [[6,7,7,8,8,9,9],[4,5,5,6,6,6,7]];
for k in [1,2] do
  xx := [t^ww[k,i]*x[i] : i in [1..7]];
  cubics := XE13Cubics(k,a,b,x);
  assert XE13Cubics(k,t^2*a,t^3*b,xx)
         eq [t^dd[k,i]*cubics[i] : i in [1..7]];
end for;

print "\nSection 5.1\n";

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);

// Matrices generating the group G = PSL_2(Z/7Z)

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;

// Functions for checking invariants and (skew) covariants.

function IsInvariant(f)
  f := PolynomialRing(K,7)!f;
  return f^M2 eq f and f^M13 eq f;
end function;

function IsCovariant(cov);
  RR := PolynomialRing(K,7);
  cov := [RR!f : f in Eltseq(cov)];
  LHS1 := [f^M2: f in cov];
  RHS1 := [&+[M2[i,j]*cov[j]: j in [1..7]]: i in [1..7]];
  LHS2 := [f^M13: f in cov];
  RHS2 := [&+[M13[i,j]*cov[j]: j in [1..7]]: i in [1..7]];
  return (LHS1 eq RHS1) and (LHS2 eq RHS2);
end function;

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]]);

function IsSkewCovariant(cov)
  RR := PolynomialRing(K,7);
  cov := [RR!f: f in Eltseq(cov)];
  LHS1 := [f^M2: f in cov];
  RHS1 := [&+[M2tilde[i,j]*cov[j]: j in [1..7]]: i in [1..7]];
  LHS2 := [f^M13: f in cov];
  RHS2 := [&+[M13tilde[i,j]*cov[j]: j in [1..7]]: i in [1..7]];
  return (LHS1 eq RHS1) and (LHS2 eq RHS2);
end function;

// The invariants of smallest degree (Q and F)

P6 := ProjectiveSpace(Rationals(),6);
R<x0,x1,x2,x3,x4,x5,x6> := CoordinateRing(P6);

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;

assert IsInvariant(Q);
assert IsInvariant(F);

print "Theorem 5.2";

cubics := [ R!f : f in cubics0 ];
H := func<f|Matrix(7,7,[Deriv(Deriv(f,i),j): i,j in [1..7]])>;
HQinv := H(Q)^(-1);
brac := func<f,g|Trace(H(f)*HQinv*H(g)*HQinv)>;

// Checking Theorem 5.2(i)

mons3 := MD(R,3);
cc := [brac(m,F - 3*Q^2) : m in mons3];
mat := Matrix(#mons3,#mons3,[MC(c,m): m in mons3,c in cc]);
km := KernelMatrix(mat);
assert Nrows(km) eq 7;
cubics1 := [&+[km[i,j]*mons3[j]: j in [1..#mons3]]: i in [1..7]];
assert HaveSameSpan(cubics,cubics1,3);

mons4 := MD(R,4);
cc := [[brac(f,m): f in cubics]:m in mons4];
mat := Matrix(#mons4,7*#mons3,[[MC(cc[i,j],m)
             : m in mons3,j in [1..7]]: i in [1..#mons4]]);
assert Nrows(KernelMatrix(mat)) eq 1;

// We already checked Theorem 5.2(ii) in the part of this
// file relating to Theorem 2.1.

// The z-curve

R6<a1,a2,a3,a4,a5,a6> := PolynomialRing(Rationals(),6);
aa := [0,a1,a2,a3,a4,a5,a6,-a6,-a5,-a4,-a3,-a2,-a1];
a := func<i|aa[(i mod 13) + 1]>;
M := Matrix(R6,13,13,[a(i-j)*a(i+j): i,j in [0..12]]);
C := Curve(Proj(R6),Pfaffians(M,4));
assert Degree(C) eq 35;

// The A-curve

printf "The scheme defined by the cubics has ";
CC := Scheme(P6,cubics);
printf "dimension %o and degree %o\n",Dimension(CC),Degree(CC);
assert Curve(P6,[x0,x2,x3,x5,x6]) subset CC;

printf "The scheme defined by the cubics and F + Q^2 has ";
CC := Scheme(P6,cubics cat [F + Q^2]);
printf "dimension %o and degree %o\n",Dimension(CC),Degree(CC);

printf "The scheme defined by the cubics, Q and F, has ";
XX := Scheme(P6,cubics cat [Q,F]);
printf "dimension %o and degree %o\n",Dimension(XX),Degree(XX);
assert IsReduced(XX);
P1 := XX![0,1,0,0,0,0,0];

// The map from the z-curve to the A-curve

ff := [1] cat [a(2*i)/a(i): i in [1,2,4,8,3,6]];
ff := [(a1*a2*a3*a4*a5*a6)*f : f in ff];
psi := map<C->CC|ff>;

print "Lemma 5.3";

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 #quartics eq 42;

P<t> := PowerSeriesRing(Integers());

X := Curve(Proj(R),cubics cat quartics);
assert Degree(X) eq 42;
mb := MinimalBasis(X);
assert {* Degree(f): f in mb *} eq {* 3^^7, 4^^42 *};
print "Hilbert series (for S) =",P!HilbertSeries(Ideal(X)) + O(t^6);
// 1 + 7*t + 28*t^2 + 77*t^3 + 119*t^4 + 161*t^5 + O(t^6)
      
XX := Scheme(Proj(R),[Q] cat cubics cat quartics);
assert Degree(XX) eq 84;
mb := MinimalBasis(XX);
assert {* Degree(f): f in mb *} eq {* 2, 3^^7, 4^^49 *};
print "Hilbert series (for T) =",P!HilbertSeries(Ideal(XX)) + O(t^6);
// 1 + 7*t + 27*t^2 + 70*t^3 + 84*t^4 + 84*t^5 + O(t^6)

// Functions for computing new invariants/covariants from old

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>;
function MyEval(cov0,cov1,cov2,m)
  P<t> := PolynomialRing(R);
  vec := Vector(P,7,[cov1[i] + t*cov2[i] : i in [1..7]]);
  cc := Compose(cov0,vec);
  return Vector(R,7,[Coefficient(cc[j],m): j in [1..7]]);
end function;

print "Lemma 5.7";

mons3 := MD(R,3);
EE := [Matrix(7,7,[<i,j,R.k>,<j,i,-R.k>]): i,j,k in [1..7]| i lt j];
dF := Vector([Deriv(F,i) : i in [1..7]]);
mat := Matrix(#EE,7*#mons3,[[MC(f,m): m in mons3, f in ff]
     where ff is [brac(vv[i],dF[i]) : i in [1..7]]
     where vv is dF*HQinv*E : E in EE]);
km := KernelMatrix(mat);
assert Nrows(km) eq 1;
N := -&+[km[1,i]*EE[i] : i in [1..#EE]];
assert N eq 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 ]
]);
vv := dF*HQinv*N;
assert forall{i : i in [1..7] | brac(vv[i],dF[i]) eq 0};

// Some functions for converting between Phi and N

function wedge1(terms)
  ii := [[i : i in [1..7] | Degree(t,i) gt 0][1]: t in terms];
  eps := 1;
  for ctr in [1..#terms] do
    for u in [1..#terms-1] do
      if ii[u] gt ii[u+1] then
        temp := ii[u]; ii[u] := ii[u+1]; ii[u+1] := temp;
        eps *:= -1;
      end if;
    end for;
  end for;
  return eps*&*terms;
end function;

function wedge(forms)
  ans := 0;
  for t1 in Terms(forms[1]) do
    for t2 in Terms(forms[2]) do
      for t3 in Terms(forms[3]) do
        ans +:= wedge1([t1,t2,t3]);
      end for;
    end for;
  end for;
  return ans;
end function;

function MatrixToForm(N)
  R := BaseRing(N);
  return (1/6)*&+[wedge([R.i,R.j,N[i,j]]): i,j in [1..7]];
end function;

function FormToMatrix(Phi)
  R := Parent(Phi);
  mat := H(Phi);
  ans := 0*mat;
  for i,j,k in [1..7] do
    xx := [i,j,k];
    Sort(~xx,~perm);
    ans[i,j] +:= Sign(perm)*MC(mat[i,j],R.k)*R.k;
  end for;
  return ans;
end function;

Phi := MatrixToForm(N);
print "Phi =",Phi;
assert N eq FormToMatrix(Phi);

print "Theorem 5.8";

v3 := nablaQ(F);
// The minus sign in the next line compensates for the fact
// the paper uses column vectors and Magma uses row vectors.
v4 := v3*(-N)*H(Q)^(-1);
v9 := Compose(v3,v3);
c6 := Dot(v4,v9);

printf "Some points on the A-curve -- above j = ";

printf "0, ";
K3<w> := CyclotomicField(3);
P<x> := PolynomialRing(K3);
L<a> := NumberField(x^2 - (-1 + 3*w));
// Fixed by M6^2 (of order 3) and so lies above j = 0
pt3 := [-2,w + a,w - a,w + a,w - a,w + a,w - a];
check := X(L)!pt3;
assert Evaluate(Q,pt3) eq -4*a^2;
assert Evaluate(c6,pt3) eq -2^16*(w - w^2)^3*a^13;
assert Evaluate(c6,pt3)^2 - 1728*Evaluate(Q,pt3)^13 eq 0;

printf "1728, ";
K4<i> := QuadraticField(-1);
P<x> := PolynomialRing(K4);
L<b> := NumberField(x^3 - (i + 1)*x^2 - x + i);
s := hom<L -> L|b^2 - b>;
// Fixed by M6^3 (of order 2) and so lies above j = 1728
pt2 := [1,b,s(b),s(s(b)),b,s(b),s(s(b))];
check := X(L)!pt2;
assert Evaluate(c6,pt2) eq 0;

print "infinity";
ptinfty := [0,1,0,0,0,0,0];
check := X!ptinfty;
assert Evaluate(Q,ptinfty) eq 0; // therefore a cusp
assert Evaluate(c6,ptinfty) eq -1;

print "\nSections 5.2 and 5.3\n";

print "Creating the covariants";

v1 := Vector([R.i : i in [1..7]]);
v3 := nablaQ(F);
v4 := v3*(-N)*H(Q)^(-1);
v9 := Compose(v3,v3);
v10 := MyEval(v4,v1,v3,3);
v12 := Compose(v3,v4);
v13 := MyEval(v4,v1,v4,3);

assert IsCovariant(v1);
assert IsCovariant(v3);
assert IsCovariant(v4);

print "Creating the skew covariants";

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;

f0 := -2*x0^3 + x0*(x1*x4 + x2*x5 + x3*x6) + x1*x3*x5 + x2*x4*x6;
f1 := x0*x1^2 + 2*x0*x3*x4 + 2*x1*x2*x6 + x2*x4^2 + x3^2*x5 + x5^2*x6;

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;

w3 := GetCovariant(f0,f1);
w4 := GetCovariant(g0,g1);
w5 := MyEval(w3,v1,v3,1);
w6 := MyEval(w4,v1,v3,1);
w7 := MyEval(w3,v1,v3,2);
w8 := MyEval(w4,v1,v3,2);
w11 := MyEval(w5,v1,v3,3);
w13 := MyEval(w5,v1,v3,4);

assert w3 eq Vector(cubics);
assert IsSkewCovariant(w3);
assert IsSkewCovariant(w4);

// g0, ... ,g6 vanish at the cusps, but not on the whole curve.
assert Eltseq(w4) subset Ideal(XX);
assert forall{g : g in Eltseq(w4) | g notin Ideal(X)};

print "Proof of Lemma 5.9 : A point over a degree 20 number field";

XX := Scheme(X,R.1-R.2);
P<x> := PolynomialRing(Rationals());
f := x^20 + 5*x^17 - 7*x^16 + 2*x^15 + 10*x^14 + x^13 + 5*x^12 +
    4*x^11 - 21*x^10 + 19*x^9 + 10*x^8 + 3*x^7 + 8*x^6 - 17*x^5 +
    5*x^3 + 2*x^2 + 1;
L<gamma> := NumberField(f);
pts := [P : P in Points(XX,L) | P[1] ne 0];
pt := [pts[1,i]/pts[1,1] : i in [1..7]];
check := X(L)!pt;

// Preparing to check the twisted versions of Q and F
// as computed in Sections 5.2 and 5.3.

RL := PolynomialRing(L,7);
S<QQ,cc6> := PolynomialRing(Rationals(),[2,13]);
RR<y1,y2,y3,y4,y5,y6,y7> := PolynomialRing(S,7);
yy := [y1,y2,y3,y4,y5,y6,y7];
Q1 := Evaluate(Q,pt);
c61 := Evaluate(c6,pt);

function GetInvariant(d,z)
  mons := MonomialsOfWeightedDegree(S,d);
  if #mons eq 0 then return 0; end if;
  mat := Matrix(#mons,20,[Eltseq(Evaluate(m,[Q1,c61])): m in mons]);
  assert Rank(mat) eq #mons;
  vec := Vector(20,Eltseq(z));
  soln := Solution(mat,vec);
  poly := &+[(Rationals()!soln[i])*mons[i] : i in [1..#mons]];
  return poly;
end function;

covs := [[v1,v3,v4,v9,v10,v12,v13],[w4,w5,w6,w7,w8,w11,w13]];
degs := [[1,3,4,9,10,12,13],[4,5,6,7,8,11,13]];

function TwistedForm(k,form,mat)
  form1 := (RL!form)^mat;
  dd := Degree(form);
  mons := MD(RR,dd);
  form2 := 0;
  for mon in mons do
    d := &+[degs[k,i]*Degree(mon,RR.i): i in [1..7]];
    c := GetInvariant(d,MC(form1,RL!mon));
    form2 +:= c*mon;
    // printf ".";
  end for;
  print "";
  return form2;
end function;

// Checking the twisted versions of Q and F as computed
// in Sections 5.2 and 5.3.

for k in [1,2] do
  casename := k eq 1 select "direct" else "skew";
  printf "Twisting the invariants in the %o case\n",casename;
  mat := Matrix(L,7,7,[Evaluate(c[i],pt): c in covs[k],i in [1..7]]);
  det := Determinant(mat);
  if k eq 1 then  // Lemma 5.9
    assert GetInvariant(52,det) eq (cc6^2 - 1728*QQ^13)^2;
  else            // Lemma 5.13
    assert GetInvariant(54,det) eq 2*QQ*(cc6^2 - 1728*QQ^13)^2;
  end if; 
  Q1,F1 := XE13RawTwistedInvariants(k,QQ,cc6,yy);
  printf "Twisting Q ";
  assert TwistedForm(k,Q,mat) eq Q1;
  printf "Twisting F ";
  assert TwistedForm(k,F,mat) eq F1;
end for;

// A procedure to check the 7 cubic forms with coefficients
// in Z[a,b] are independent. The coefficients of the cubic
// forms may be arranged in a 7 by 84 matrix. We check that
// the 7 by 7 minors generate the unit ideal in Z[a,b].
// Doing this by choosing the minors at random can take
// a little time, so we have coded in which minors to use.

procedure CheckIndependent(forms)
  R := Universe(forms);
  d := Degree(forms[1]);
  assert forall{f : f in forms | Degree(f) eq d};
  mons := MD(R,d);
  mat := Matrix(#forms,#mons,[MC(f,m): m in mons,f in forms]);
  P<A,B> := PolynomialRing(Integers(),2);
  mat := Matrix(P,#forms,#mons,[Numerator(mat[i,j])
                : j in [1..Ncols(mat)], i in [1..Nrows(mat)]]);
  n := #forms;
  ii := [];
  if [d,n] eq [3,7] then
    ii := [
      [ 14, 31, 50, 53, 59, 67, 81 ],
      [ 17, 18, 26, 32, 41, 45, 79 ],
      [ 17, 29, 31, 32, 36, 53, 84 ],
      [  3, 17, 53, 61, 71, 75, 77 ],
      [ 28, 49, 50, 54, 74, 80, 84 ],
      [ 28, 53, 54, 62, 74, 83, 84 ]];
  end if;
  g := 0;
  ctr := 0;
  repeat
    ctr +:= 1;
    if ctr le #ii then
      jj := ii[ctr];
    else
      jj := Sort([Random([1..#mons]): i in [1..n]]);
    end if;
    submat := Matrix(P,n,n,[mat[i,j] : j in jj,i in [1..n]]);
    det := Determinant(submat);
    if det ne 0 then
      g := GCD(g,Numerator(det));
    end if;
  until g eq 1;
end procedure;

// The changes of coordinates used to simplify the twisted
// forms of Q and F. (This is needed to obtain correct
// formulae when j(E) = 0.)

K<a,b> := FunctionField(Rationals(),2);
c6 := -864*b;
P<X> := PolynomialRing(K);
L<Q> := quo<P|X^13 - 16*(4*a^3 + 27*b^2)>;
RL<x1,x2,x3,x4,x5,x6,x7> := PolynomialRing(L,7);

yy := [[ 16*Q^6*(12*a^2*(x1 + 2*x2) + 18*a*b*x3
    + 14*(4*a^3 + 27*b^2)*x5 + 81*b^2*x6),
  2*Q^5*(96*a^2*x2 - 144*a*b*(x3 - 3*x4) + (48*a^3 - 324*b^2)*x5
    + (52*a^3 - 297*b^2)*x6  - 864*a^2*b*x7),
  Q^11*(2*a*(x3 - 56*x4) + 3*b*(44*x5 + 3*x6) + 176*a^2*x7),
  -4*Q^2*(4*a^3 + 27*b^2)*(2*x5 - x6),
  2*Q^8*(2*a*x4 - 3*b*x5 - 4*a^2*x7),
  Q^7*(2*a*(x3 - 2*x4) + 3*b*(2*x5 + 3*x6) + 8*a^2*x7),
  2*(4*a^3 + 27*b^2)*x6 ],[
  6*Q^11*(36*b*x1 + 90*b*x2 - 2*a^2*(5*x3 + 8*x4 + 2*x5)
    - 9*a*b*x6 + 42*a*b*x7),
  4*Q^4*(4*a^3 + 27*b^2)*(1347*x1 - 936*x2 - 317*a*x6 - 185*a*x7),
  3*Q^10*(36*b*x1 - 18*b*x2 + 2*a^2*(x3 - 8*x4 - 14*x5)
     + 9*a*b*x6 + 6*a*b*x7),
  24*Q^3*(4*a^3 + 27*b^2)*(21*x1 - 36*x2 - a*x6 + 9*a*x7),
  3*Q^9*(18*b*x2 - 2*a^2*(x3 + 2*x5) + 3*a*b*x6 + 6*a*b*x7),
  8*Q*(4*a^3+27*b^2)*(6*x1 - 6*x2 - a*x6),
  -4*(4*a^3+27*b^2)*(3*x1 - a*x6 - a*x7) ]];

det := [-2^23*3^3*a^8*Q^39*(4*a^3 + 27*b^2)^2,
         2^26*3^8*a^8*Q^38*(4*a^3 + 27*b^2)^4];
la := [2^16*3^2*a^4*(4*a^3 + 27*b^2),
         2^18*3^4*a^4*(4*a^3 + 27*b^2)^2];
mu := [2^32*3^4*a^8*(4*a^3 + 27*b^2)^2,
         2^32*3^8*a^8*(4*a^3 + 27*b^2)^3];

RK<x1,x2,x3,x4,x5,x6,x7> := PolynomialRing(K,7);
xx := [RK.i : i in [1..7]];

for k in [1,2] do

  casename := k eq 1 select "direct" else "skew";
  printf "Simplifying the twisted invariants in the %o case\n",casename;
  M := Matrix(7,7,[MC(f,RL.i): i in [1..7],f in yy[k]]);
  assert Determinant(M) eq det[k];
  Q,F := XE13RawTwistedInvariants(k,L.1,c6,yy[k]);
  Q := (RK!Q)/la[k];
  F := (RK!F)/mu[k];
  Q1,F1 := XE13TwistedInvariants(k,a,b,xx);
  assert Q eq Q1 and F eq F1;

// Checking the cubics in congr13-formulae.m span the spaces of
// cubics described in Theorem 5.10 and 5.14.

  cubics := XE13Cubics(k,a,b,xx);
  assert #cubics eq 7;
  HQinv := H(Q)^(-1);
  brac := func<f,g|Trace(H(f)*HQinv*H(g)*HQinv)>;
  rho := k eq 1 select 1 else 16*(4*a^3 + 27*b^2);
  assert forall{f : f in cubics | brac(f,F - 3*rho*Q^2) eq 0};
  CheckIndependent(cubics);

// Checking the alternating matrices N in congr13-formulae.m
// satisfy the analogue of Lemma 5.7 where Q and F are
// replaced by their twisted versions.

  N := XE13AlternatingMatrix(k,a,b,xx);
  dF := Vector([Deriv(F,i) : i in [1..7]]);
  ww := dF*HQinv*N;
  assert forall{i : i in [1..7] | brac(ww[i],dF[i]) eq 0};
  Phi := MatrixToForm(N);
  print "Phi =",Phi;
  assert N eq FormToMatrix(Phi);
  assert GCD(Pfaffians(N,6)) eq Q;

// Checking Remark 5.11

  if k eq 1 then
    subst := [
      -a*x1 - a^2*x2 + 9*b*x3 + 9*b*x4 - (3/2)*b*x5
         - 3*a*x6 + 2*a^2*x7,
      (1/2)*a*x1 + (1/3)*a^2*x2 - (9/2)*b*x3 - 3*b*x4
         + (3/4)*b*x5 + a*x6 - a^2*x7,
      3*b*x2 + 2*a*x4, 2*a*x3 - 3*b*x7, -(1/2)*x1 - a*x7,
      -(1/3)*a*x2 + x6, (1/2)*x3 + (1/12)*x5 ];
    Q0 := Evaluate(Q,subst)/(4*a^3 + 27*b^2);
    N0 := Evaluate(N,subst);
    cob := Matrix(7,7,[MC(f,RK.i): i in [1..7],f in subst]);
    N0 := (1/(4*a^3 + 27*b^2))*Transpose(cob)*N0*cob;
    Phi0 := MatrixToForm(N0);
    assert Q0 eq x1*x7 + x2*x6 + x3*x5 + x4^2;
    assert Phi0 eq x1*x2*x3 + x1*x4*x7 + x2*x4*x6 + x3*x4*x5 + x5*x6*x7;
    assert N0 eq FormToMatrix(Phi0);
 end if;
end for;

// Function creating X_E(13,k) as a scheme in Magma
// using the formulae in Theorems 5.10 and 5.14.

function XE13twist(k,ab)
  assert k in [1,2];
  K0 := Universe(ab);
  if K0 cmpeq Integers() then K0 := Rationals(); end if;
  R := PolynomialRing(K0,7);
  xx := [R.i : i in [1..7]];
  a,b := Explode(ab);
  rho := k eq 1 select 1 else 16*(4*a^3 + 27*b^2);
  Q,F := XE13TwistedInvariants(k,a,b,xx);
  cubics := XE13Cubics(k,a,b,xx);
  return Scheme(Proj(R),cubics cat [F + rho*Q^2]);
end function;

// Function evaluating the j-map X_E(13,k) -> P1

function jmapXE13(k,ab,pt);
  assert k in [1,2];
  K0 := Universe(ab);
  if K0 cmpeq Integers() then K0 := Rationals(); end if;
  R := PolynomialRing(K0,7);
  xx := [R.i : i in [1..7]];
  a,b := Explode(ab);
  Q,F := XE13TwistedInvariants(k,a,b,xx);
  N := XE13AlternatingMatrix(k,a,b,xx);
  HQ := H(Q);
  HQinv := HQ^(-1);
  v3 := Vector([Deriv(F,i) : i in [1..7]])*HQinv;  
  v4 := -v3*N*HQinv;
  pt := Eltseq(pt);
  v3pt := Evaluate(v3,pt);
  v4pt := Evaluate(v4,pt)[1];
  v9pt := Evaluate(v3,Eltseq(v3pt))[1];
  QQ := Evaluate(Q,pt);
  cc6 := &+[HQ[i,j]*v4pt[i]*v9pt[j] : i,j in [1..7]];
  if k eq 1 then
    return 1728 - (4*a^3 + 27*b^2)*cc6^2/QQ^13;
  else
    return 1728 - cc6^2/(2^40*(4*a^3 + 27*b^2)^10*QQ^13);
  end if;
end function;

print "\nSection 5.4\n";

// Creating the 2-isogenous elliptic curves E and F

K<a,b> := FunctionField(Rationals(),2);
P<X> := PolynomialRing(K);
L<theta> := quo<P|X^3 + a*X + b>;
A := -15*theta^2 - 4*a;
B := 14*a*theta + 22*b;
E := EllipticCurve([L|a,b]);
F := EllipticCurve([A,B]);
c := 3*theta;
d := 3*theta^2 + a;
assert IsIsomorphic(E,EllipticCurve([0,c,0,d,0]));
assert IsIsomorphic(F,EllipticCurve([0,-2*c,0,c^2-4*d,0]));

// Checking that X_E(13,1) and X_E(F,2) are isomorphic

RL<x1,x2,x3,x4,x5,x6,x7> := PolynomialRing(L,7);
xx := [RL.i : i in [1..7]];
cubics1 := XE13Cubics(1,a,b,xx);
cubics2 := XE13Cubics(2,A,B,xx);

subst := [
  (14*theta^2 + 8*a)*x1 + 32*theta^2*x2 - (3*a*theta - 25*b)*x3
    + (38*a*theta - 18*b)*x4 + (-18*a*theta^2 - 30*b*theta - 16*a^2)*x5
    - (11*a*theta^2 - 3*b*theta + 24*a^2)*x6
    + (144*b*theta^2 - 44*a^2*theta + 132*a*b)*x7,
  (14*theta^2 + 8*a)*x2 + (3*a*theta + 7*b)*x3 + (11*a*theta + 15*b)*x4
    - (17*a*theta^2 - 9*b*theta + 8*a^2)*x5
    - (5*a*theta^2 - 21*b*theta + 8*a^2)*x6
    + (144*b*theta^2 + 10*a^2*theta + 66*a*b)*x7,
  8*theta*x1 + 8*theta*x2 + 8*theta^2*x3 + (30*theta^2 + 24*a)*x4
    + (10*a*theta - 6*b)*x5 - (8*a*theta + 24*b)*x6
    + (12*a*theta^2 - 108*b*theta)*x7,
  4*theta*x1 + 8*theta*x2 + (5*theta^2 + 4*a)*x3 - (5*a*theta - 3*b)*x6,
  -4*theta*x1 - 12*theta*x2 + (3*theta^2 - 4*a)*x4 + (a*theta + 9*b)*x5
    + (6*a*theta^2 - 18*b*theta + 16*a^2)*x7,
  -4*x1 - 12*x2 + 6*theta*x4 - (6*theta^2 + 4*a)*x5 - 12*a*theta*x7,
  2*x1 + 4*x2 + theta*x3 + (3*theta^2 + 2*a)*x6
];

mat := Matrix(7,7,[MC(subst[i],RL.j): i,j in [1..7]]);
assert Determinant(mat) eq -2^10*3^2*d^3*(c^2 - 4*d)^5;
cubics2a := [Evaluate(f,subst) : f in cubics2];
assert HaveSameSpan(cubics1,cubics2a,3);
mons := MD(RL,3);
mat1 := Matrix(7,#mons,[MC(f,mon): mon in mons,f in cubics1]);
mat2 := Matrix(7,#mons,[MC(f,mon): mon in mons,f in cubics2a]);
cob := Solution(mat1,mat2);
assert Determinant(cob) eq 2^45*d^6*(c^2 - 4*d)^18;

// Checking the theorems hold for y^2 = x^3 + x and y^2 = x^3 + 1

E10 := EllipticCurve([1,0]);
F := EllipticCurve([-44,112]);
flag,isog := IsIsogenous(E10,F);
assert flag and Degree(isog) eq 4;

E01 := EllipticCurve([0,1]);
F := EllipticCurve([-15,22]);
flag,isog := IsIsogenous(E01,F);
assert flag and Degree(isog) eq 2;

// Checking the final paragraph of Section 5.4

E := EllipticCurve([0,2]);
F := EllipticCurve([-120,506]);
flag,isog := IsIsogenous(E,F);
assert flag and Degree(isog) eq 3;

for p in [17,19,23,37] do
  X1 := XE13twist(1,[GF(p)|0,2]);
  X2 := XE13twist(1,[GF(p)|0,1/2]);
  X3 := XE13twist(1,[GF(p)|-120,506]);
  print p,[#Points(X): X in [X1,X2,X3]];
end for;

print "\nSection 8\n";

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

print "Example 8.1";
a := -4;
b := -3;
E1 := EllipticCurve([a,b]);
assert CremonaReference(E1) eq "52a2";
C := XE13twist(1,[a,b]);
assert Degree(C) eq 42;
P1 := C![1,0,0,0,0,0,0];
P2 := C![-30,23,-72,-16,0,16,1];
E2 := EllipticCurve("988b1");
assert jInvariant(E2) eq -2^8*3^3*151^3*2399^3/(13*19^13);
assert jmapXE13(1,[a,b],P1) eq jInvariant(E1);
assert jmapXE13(1,[a,b],P2) eq jInvariant(E2);

print "Example 8.2";
a := 1;
b := -10;
E0 := EllipticCurve([a,b]);
assert CremonaReference(E0) eq "52a1";
C := XE13twist(2,[a,b]);
assert Degree(C) eq 42;
P1 := C![9,0,4,2,-2,2,-1];
P2 := C![134,-45,134,44,5,-18,4];
E1 := EllipticCurve("52a2");
E2 := EllipticCurve("988b1");
assert jInvariant(E1) eq 2^14*3^3/13;
assert jInvariant(E2) eq -2^8*3^3*151^3*2399^3/(13*19^13);
assert jmapXE13(2,[a,b],P1) eq jInvariant(E1);
assert jmapXE13(2,[a,b],P2) eq jInvariant(E2);

print "Example 8.5";
K<t> := FunctionField(Rationals());
p := t^2 - 3*t - 1;
q := t^3 - 2*t^2 - 3*t - 8;
r := t^3 + t^2 + 2*t - 1;
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;
a := -3*p^2*q*f1;
b := 2*p^2*q^2*f2;
EE := EllipticCurve([a,b]);
C := XE13twist(1,[a,b]);
pt := [ 2592*t*(t + 1)*(t^2 + 2*t + 3)*p^2*q*r^2,
   -432*t^2*(2*t^2 + 3*t + 7)*p^2*q*r^2,
   -72*(t^3 + 2*t^2 + 4*t - 1)*p^2*q*r,
   72*p^2*q*r, 
   -6*(t^6 - 3*t^4 - 11*t^3 - 14*t^2 - 5*t - 4)*p,
   -24*(t - 1)*p*r,
   t - 1 ];
check := C!pt;
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;
d := t^4*(t + 2)^3*(t^4 + 4*t^3 + 9*t^2 + 11*t + 8)*(t^6 + 4*t^5 +
  9*t^4 + 8*t^3 + 2*t^2 - 9*t - 6)^13;
j1 := p*g1^3/(r*d);
assert j1 eq 1728 + q*g2^2/(r*d);
assert jmapXE13(1,[a,b],pt) eq j1;
EE1 := EllipticCurve([-3*p^3*q^3*g1,2*p^4*q^5*g2]);
assert jInvariant(EE1) eq j1;

Table2 := [
  <1, -6, "11a3", "11a2", 25>,
  <-1, -1, "768h1", "768h4", 10>,
  <4, 3, "13688b1", "8363368*", 0>,
  <2, 2, "27930s1", "27930r1", 0>,
  <-4, 1, "80408l1", "8282024*", 0>,
  <-1/2, -6, "83030b1", "913330*", 0>, 
  <-1/3, 14, "271545f1", "589524195*", 0>,
  <1/2, 6, "5429670*", "320350530*", 0>,
  <1/4, -1, "7707798*", "27925352154*", 0>,
  <-3, 6, "15211515*", "1566786045*", 0>,
  <8/5, -5, "46427580*", "5448415795740*", 0>,
  <3, -3, "48963840*", "42941287680*", 0>,
  <5, 6, "147656145*", "1624217595*", 0>,
  <7/2, 2, "192105606*", "23030964786522*", 0>,
  <5/2, -2, "774703710*", "6034167197190*", 0>,
  <7, -1, "1040014080*", "24181367374080*", 0> ];

for dd in Table2 do
  t0,d,r1,r2,ii := Explode(dd);
  E := EllipticCurve([Evaluate(a,t0): a in Eltseq(EE)]);
  E1 := EllipticCurve([Evaluate(a,t0): a in Eltseq(EE1)]);
  E := MinimalModel(QuadraticTwist(E,d));
  E1 := MinimalModel(QuadraticTwist(E1,d));
  score := #[x : x in Eltseq(aa(E) - aa(E1)) | x ne 0];
  assert score le 2;
  assert CR(E) eq r1;
  assert CR(E1) eq r2;
  flag,isog := IsIsogenous(E,E1);
  assert ii eq (flag select Degree(isog) else 0);
end for;

print "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;
a := -3*p*q^2*f1;
b := 2*p*q^2*f2;
EE := EllipticCurve([a,b]);
r1 := 29*t^7 + 15*t^6 + 7*t^5 - 32*t^4 + 89*t^3 + 73*t^2 + 83*t + 24;
r2 := 55*t^7 + 117*t^6 + 269*t^5 + 356*t^4 + 211*t^3
      + 179*t^2 + 25*t + 36;
r3 := t^6 - 14*t^5 - 43*t^4 - 85*t^3 - 106*t^2 - 53*t - 36;
r4 := 13*t^5 + 34*t^4 + 17*t^3 + 11*t^2 - 22*t - 5;
pt := [ -2*p*q*r1, p*q*r2, 2*r3, -24*(t + 1)^2*(t^2 + 1)^2,
        t*r4, 2*t, 0 ];
C := XE13twist(2,[a,b]);
check := C!pt;
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;
d := (t+1)*(t^2+1)*(2*t^2 + t + 1)^2*(2*t^3 + 2*t^2 + 3*t + 1)^13;
j1 := p*q*g1^3/d;
assert j1 eq 1728 - g2^2/d;
assert jmapXE13(2,[a,b],pt) eq j1;
EE1 := EllipticCurve([3*p*q^3*g1,2*p*q^4*g2]);
assert jInvariant(EE1) eq j1;

Table3 := [
  <0, 3, "121c1", "121c2", 11>,
  <1, 3, "162c1", "162c4", 21>,
  <-1/3, -10, "1225h1", "1225h2", 37>,
  <-3, -21, "1960i1", "21560l1", 0>,
  <-2, -1, "14175k1", "184275o1", 0>,
  <1/3, -26, "23660f1", "733460*", 0>,
  <-3/4, -39, "92950q1", "2881450*", 0>,
  <-1/2, 22, "98010s1", "98010t1", 0>,
  <3, 78, "185900a1", "7621900*", 0>,
  <-7, 645, "255162e1", "4848078*", 0>,
  <1/2, -30, "1242150*", "1242150*", 0>,
  <1/7, 210, "1695978*", "429082434*", 0>,
  <-3/2, -6, "2141594*", "49256662*", 0>,
  <-3/5, -3, "2147950*", "2147950*", 0>,
  <-5, 6, "4746924*", "507920868*", 0>,
  <1/5, 30, "7495800*", "397277400*", 0> ];

for dd in Table3 do
  t0,d,r1,r2,ii := Explode(dd);
  E := EllipticCurve([Evaluate(a,t0): a in Eltseq(EE)]);
  E1 := EllipticCurve([Evaluate(a,t0): a in Eltseq(EE1)]);
  E := MinimalModel(QuadraticTwist(E,d));
  E1 := MinimalModel(QuadraticTwist(E1,d));
  score := #[x : x in Eltseq(aa(E) - aa(E1)) | x ne 0];
  assert score le 2;
  assert CR(E) eq r1;
  assert CR(E1) eq r2;
  flag,isog := IsIsogenous(E,E1);
  assert ii eq (flag select Degree(isog) else 0);
end for;

print "";
