
// Title : On pairs of 17-congruent elliptic curves
// Author : T. Fisher
// Version : June 2021

// This is a MAGMA file checking some of the calculations in
// the paper, in roughly the same order as the paper.

// We start by loading in some formulae from the file

load "congr17-formulae.m";

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

// Some abbreviations

MC := MonomialCoefficient;
Deriv := Derivative;
Coeff := Coefficient;
EC := EllipticCurve;
N := Conductor;
a := TraceOfFrobenius;
aa := TracesOfFrobenius;
IP := InnerProduct;
Symm := Symmetrization;
QQ := Rationals();
QQT<T> := FunctionField(QQ);

printf "\nSection 1 : Introduction\n\n";

// The new pair of 17-congruent elliptic curves

E1 := EC([1,-1,0,-128973503459,17827877649739965]);
E2 := EC([1,-1,0,-184201215542543714,-34187608332483214491862380]);

assert N(E1) eq 2*5*13*59*191^2;
assert N(E2) eq 2*5*11*13*59*191^2;

// The old pair of 17-congruent elliptic curves

E1dash := EC([1,0,0,-8,27]);
E2dash := EC([1,0,0,8124402,-11887136703]);

assert N(E1dash) eq 3*5^2*7^2;
assert N(E2dash) eq 3*5^2*7^2*13;

E := EC([ T^2 - T - 2, -1, T^3, 0, 0 ]);
assert Regulator([E|[0,0],[-T,-T]]) eq 13/36;

printf "Section 2 : Verification via modularity\n\n";

mu := func<M|Integers()!(M*&*[1+1/l : l in PrimeDivisors(M)])>;

assert aa(E1,50) eq [-1,0,1,-1,-5,1,2,4,-1,-3,-2,-11,5,-4,-9];
assert aa(E2,50) eq [-1,0,1,-1,1,1,2,4,-1,-3,-2,6,-12,-4,-9];
sturmbound := mu(M)/6 where M is N(E2);
assert sturmbound eq 1.0328006744384765625*2^30;
p := 11;
repeat
  p := NextPrime(p);
  assert GF(17)!(a(E1,p) - a(E2,p)) eq 0;
until p gt 10^5; // sturmbound;
assert Discriminant(E1) eq 2^3*5^3*13*59^2*191^3;
assert Discriminant(E2) eq -2^14*5^11*11^17*13*59*191^9;

assert aa(E1dash,50) eq [-1,1,0,0,0,-3,2,-1,-2,-8,8,-7,0,8,-10];
assert aa(E2dash,50) eq [-1,1,0,0,0,1,2,-1,-2,9,-9,10,0,8,7];
sturmbound := mu(M)/6 where M is N(E2dash);
assert sturmbound eq 15680;
p := 13;
repeat
  p := NextPrime(p);
  assert GF(17)!(a(E1dash,p) - a(E2dash,p)) eq 0;
until p gt sturmbound;
assert Discriminant(E1dash) eq -3^5*5^2*7^2;
assert Discriminant(E2dash) eq -3^2*5^2*7^2*13^17;

printf "Section 3 : Verification via genus 2 Jacobians\n\n";
IndentPush();

Cprint<i> := ComplexField(10);
prec := 300;
tau1 := w[2]/w[1] where w is Periods(E1:Precision := prec);
tau2 := w[2]/w[1] where w is Periods(E2:Precision := prec);

assert Sprint(Cprint!tau1) eq "0.1142862335*i";
assert Sprint(Cprint!tau2) eq "0.5000000000 + 1.415897663*i";

tau2dash := (64*tau2 - 15)/(-17*tau2 + 4);
tau := Matrix(2,2,[17*tau1,tau1,tau1,(tau1 + tau2dash)/17]);

print "Computing the Rosenhain invariants";

time rr := RosenhainInvariants(tau);
Cx<x> := PolynomialRing(ComplexField(prec));
f0 := x*(x-1)*&*{x-a : a in rr};
IC := IgusaClebschInvariants(f0);
IC := [IC[i]/IC[1]^[1,2,3,5][i] : i in [1..4]];
IC0 := [BestApproximation(Re(r),10^150) : r in IC];

print "Computing the genus 2 curve";

time C0 := HyperellipticCurveFromIgusaClebsch(IC0);
time C0 := ReducedWamelenModel(C0);

QQx<x> := PolynomialRing(QQ);
f1,f2,g1,g2,h1,h2 := Genus2Polynomials(QQx);
assert f1 eq 196081931*x^3 + 1143338037*x^2 - 801791940*x + 135616700;
assert f2 eq -25996*x^3 + 1698260*x^2 - 6845267*x + 3822078;
C := HyperellipticCurve(f1*f2);
IC := IgusaClebschInvariants(C);
IC := [IC[i]/IC[1]^[1,2,3,5][i] : i in [1..4]];
assert IC eq IC0;
assert IsQuadraticTwist(C,C0);

for p in [240..1000] do
  if IsPrime(p) then
    Fpx := PolynomialRing(GF(p));
    Cp := HyperellipticCurve(Fpx!(f1*f2));
    E1p := EC([GF(p)!a : a in Eltseq(E1)]);
    E2p := EC([GF(p)!a : a in Eltseq(E2)]);
    assert #Jacobian(Cp) eq #E1p*#E2p;
  end if;
end for;

print "Checking the maps C -> E1 and C -> E2";

KC<u,v> := FunctionField(C);
xmap := [h1/(f1*g1^2),h2/(f2*g2^2)];
elldiff := [(273857*u - 336364)/v,(2758*u + 1630)/v];
for i in [1..2] do
  xi := Evaluate(xmap[i],u);
  dxi := Evaluate(Deriv(xmap[i]),u);
  eta := (dxi/elldiff[i] - xi)/2;
  psi := map<C->[E1,E2][i]|[xi,eta,1]>;
  assert Degree(psi) eq 17;
end for;

Fpx<x> := PolynomialRing(GF(101));
assert Fpx!g1 eq 25*x^7 + 56*x^6 + 31*x^5 + 99*x^4 + 100*x^3 + 42*x^2
  + 79*x + 5;
assert Fpx!g2 eq 3*x^7 + 76*x^6 + 44*x^5 + 97*x^4 + 52*x^3 + 38*x^2 +
  75*x + 2;
assert Fpx!h1 eq 16*x^17 + 6*x^16 + 57*x^15 + 54*x^14 + 94*x^13 +
  79*x^12 + 77*x^11 + 55*x^10 + 74*x^9 + 78*x^8 + 97*x^7 + 79*x^6 +
  25*x^5 + 96*x^4 + 98*x^3 + 46*x^2 + 4*x + 99;
assert Fpx!h2 eq 67*x^17 + 25*x^16 + x^15 + 22*x^14 + 84*x^13 +
  94*x^12 + 93*x^11 + 95*x^10 + 34*x^9 + 40*x^8 + 99*x^7 + 84*x^6 +
  43*x^5 + 12*x^4 + 59*x^3 + 13*x^2 + 26*x + 98;

IndentPop();

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

printf "\nSection 4 : The modular curve X(17)\n\n";
IndentPush();

K<z17> := CyclotomicField(17);
sqrt17 := &+[LegendreSymbol(a,17)*z17^a : a in [0..16]];
assert sqrt17 eq 1/z17*&*[(1 - z17^i): i in [1..8]];
assert sqrt17^2 eq 17;
xi := [z17^i + z17^(-i): i in [1..8]];
I9 := IdentityMatrix(K,9);
M2 := (-1/sqrt17)*Matrix(K,9,9,[
    1,   1,     1,     1,     1,     1,     1,     1,     1,
    2, xi[3], xi[8], xi[7], xi[4], xi[5], xi[2], xi[6], xi[1],
    2, xi[8], xi[7], xi[4], xi[5], xi[2], xi[6], xi[1], xi[3],
    2, xi[7], xi[4], xi[5], xi[2], xi[6], xi[1], xi[3], xi[8],
    2, xi[4], xi[5], xi[2], xi[6], xi[1], xi[3], xi[8], xi[7],
    2, xi[5], xi[2], xi[6], xi[1], xi[3], xi[8], xi[7], xi[4],
    2, xi[2], xi[6], xi[1], xi[3], xi[8], xi[7], xi[4], xi[5],
    2, xi[6], xi[1], xi[3], xi[8], xi[7], xi[4], xi[5], xi[2],
    2, xi[1], xi[3], xi[8], xi[7], xi[4], xi[5], xi[2], xi[6]
]);
M8 := Matrix(9,9,[-I9[j] : j in [1,9,2,3,4,5,6,7,8]]);
M17 := DiagonalMatrix(K,[z17^i : i in [0,1,9,13,15,16,8,4,2]]);
assert forall{M : M in [M2,M8,M17] | Determinant(M) eq 1};

G := MatrixGroup<9,K|M2,M17>;
assert #G eq 2448;
assert IsIsomorphic(G,PSL(2,17));
assert M8 in G;
M3 := M17^11*M2*M17^10*M8^2;
assert Order(M3) eq 3;
sigma := hom<K->K|z17^3>;
M2tilde := Matrix(9,9,[sigma(M2[i,j]): i,j in [1..9]]);
M17tilde := Matrix(9,9,[sigma(M17[i,j]): i,j in [1..9]]);

IsInvariant := func<f0| f^M2 eq f and f^M17 eq f
  where f is PolynomialRing(K,9)!f0>;
R9<x0,x1,x2,x3,x4,x5,x6,x7,x8> := PolynomialRing(QQ,9);
Q := x0^2 + x1*x5 + x2*x6 + x3*x7 + x4*x8;
D := 2*x0*(x1*x5 - x2*x6 + x3*x7 - x4*x8) - x1^2*x4 + x2^2*x5 -
  x3^2*x6 + x4^2*x7 - x5^2*x8 + x6^2*x1 - x7^2*x2 + x8^2*x3;
F := x0^4 + x0*(x1^2*x4 + x2^2*x5 + x3^2*x6 + x4^2*x7 + x5^2*x8 +
  x1*x6^2 + x2*x7^2 + x3*x8^2) + x1*x3*x5*x7 + x2*x4*x6*x8 +
  x1*x2*x5*x6 + x2*x3*x6*x7 + x3*x4*x7*x8 + x1*x4*x5*x8 + x1^2*x3*x8 +
  x1*x2^2*x4 + x2*x3^2*x5 + x3*x4^2*x6 + x4*x5^2*x7 + x5*x6^2*x8 +
  x1*x6*x7^2 + x2*x7*x8^2;
assert forall{f : f in [Q,D,F] | IsInvariant(f)};

print "Creating the z-curve";

R8<z1,z2,z3,z4,z5,z6,z7,z8> := PolynomialRing(QQ,8);
zz := [0,z1,z2,z3,z4,z5,z6,z7,z8,-z8,-z7,-z6,-z5,-z4,-z3,-z2,-z1];
z := func<i|zz[(i mod 17) + 1]>;
mat := Matrix(R8,17,17,[z(i-j)*z(i+j): i,j in [0..16]]);
eqns := Pfaffians(mat,4);
subst := [zz[(3*i) mod 17 + 1] : i in [1..8]];
MM := Matrix(QQ,8,8,[MC(subst[i],R8.j): i,j in [1..8]]);
extra_quartics := [
  z1^3*z2 + z1*z4^2*z5 + z1*z5*z6*z8 - z2*z4*z6^2,
  z1^2*z2*z8 - z1^2*z6*z7 - z1*z5^2*z6 + z2*z3*z5*z7,
  z1*z2*z6*z7 - z1*z3*z4*z8 - z3^2*z6^2 + z3*z4^2*z7,
  z1*z3*z4*z7 + z2*z3^2*z6 - z3*z6*z7^2 - z4*z5*z6*z7,
  z1*z3*z7*z8 + z1*z4*z5*z8 - z3*z4*z7^2 + z3*z5^2*z8,
  z1^2*z5*z7 - z1*z2^2*z4 - z2^2*z3*z5 + z2*z3*z4*z8 + z3^3*z7,
  z1^2*z4^2 + z1^2*z6*z8 + z1*z3*z4*z5 - z2*z3*z5*z8 - z2*z4^2*z7,
  z1^2*z5^2 - z1*z2^2*z3 + z2^2*z5*z6 - z2*z5^2*z7 + z3^3*z5 - z3*z4*z5*z6
];
eqns cat:= [f^(MM^j): j in [0..7],f in extra_quartics];
eqns := MinimalBasis(ideal<R8|eqns>);
assert {* TotalDegree(f): f in eqns *} eq {* 4^^126 *};
C := Curve(Proj(R8),eqns);
assert Degree(C) eq 84;

print "Proposition 4.1 : The curves C1 and C2";

eqns0 := [Q] cat [Deriv(F,i): i in [1..9]];
phi1 := [ 1, z2/z1, z6/z3, -z1/z8, -z3/z7, z8/z4, -z7/z5, z4/z2, -z5/z6 ];
eqns1 := eqns0 cat [x0^4 + x1*x3*x5*x7,x0^4 + x2*x4*x6*x8];
assert [Numerator(Evaluate(f,phi1)): f in eqns1] subset Ideal(C);
phi2 := [
    z1*z4 + z2*z8 + z3*z5 - z6*z7 - 2*z4*z7*(z2/z1),
    -z8^2 - 2*z5*z7*(z2/z1), z7^2 - 2*z2*z4*(z6/z3),
    -z4^2 + 2*z5*z6*(z1/z8), z5^2 + 2*z1*z2*(z3/z7),
    -z2^2 + 2*z3*z6*(z8/z4), z6^2 - 2*z1*z8*(z7/z5),
    -z1^2 - 2*z3*z7*(z4/z2), z3^2 + 2*z4*z8*(z5/z6)
];
eqns2 := eqns0 cat [D];
assert [Numerator(Evaluate(f,phi2)): f in eqns2] subset Ideal(C);

// Additional equations for C1

M8 := ChangeRing(M8,QQ);
quartics := ExtraQuarticsForACurve(R9);
quartics := [q^(M8^j): j in [0..7],q in quartics];
// I1 := ideal<R9|eqns1>;
// time ii := [i : i in [1..#quartics] | quartics[i]^3 in I1];
// I1 := ideal<R9|eqns1 cat [quartics[i] : i in ii]>;
// time assert forall{q : q in quartics | q^3 in I1};
I1 := ideal<R9|eqns1 cat quartics>;
C1 := Scheme(Proj(R9),I1);
time assert Dimension(C1) eq 1; 
assert Degree(C1) eq 96;
eqns1 := MinimalBasis(I1);
assert {* TotalDegree(f): f in eqns1 *} eq {* 2, 3^^9, 4^^117 *};
X17eqns := eqns1;

// Additional equations for C2

cubics := [
  x0^3 - x0*x1*x5 + x0*x2*x6 + x2*x7^2 + x3^2*x6,
  x0*x1*x2 + x0*x6*x8 + x1^2*x6 - x2*x5*x8 - x3*x4*x5,
  x0^2*x1 + x0*x2^2 + x0*x5*x8 + x1*x3*x7 + x2*x3^2 + x4*x5*x7
];
cubics := [c^(M8^j): j in [0..7],c in cubics];
// I2 := ideal<R9|eqns2>;
// time assert forall{f : f in cubics | f^3 in I2};
I2 := ideal<R9|eqns2 cat cubics>;
C2 := Scheme(Proj(R9),I2);
time assert Dimension(C2) eq 1;
assert Degree(C2) eq 168;
eqns2 := MinimalBasis(I2);
assert {* TotalDegree(f): f in eqns2 *} eq {* 2, 3^^28 *};

// print "Proving the decomposition C = C1 union C2";
// I := ideal<R9|eqns0>;
// time assert D^2*(x0^4 + x1*x3*x5*x7) in I;
// time assert D^2*(x0^4 + x2*x4*x6*x8) in I;

print "Lemma 4.4 : The map j : X(17) -> P^1";

H := func<f|Matrix(9,9,[Deriv(Deriv(f,i),j): i,j in [1..9]])>;
_<t> := PolynomialRing(R9);
HQ := H(Q);
HQinv := HQ^(-1);
nabla := func<f|Vector([Deriv(f,i) : i in [1..9]])*HQinv>;
v1 := Vector([R9.i : i in [1..9]]);
v2 := nabla(D);
v3 := nabla(F);
v4 := Vector([Evaluate(f,Eltseq(v2)): f in Eltseq(v2)]);
v6 := Vector([Evaluate(f,Eltseq(v2)): f in Eltseq(v3)]);
cc4 := &+[HQ[i,j]*v4[i]*v6[j] : i,j in [1..9]];
assert Evaluate(cc4,[0,1,0,0,0,0,0,0,0]) eq 0;
Qi<i> := QuadraticField(-1);
P<X> := PolynomialRing(Qi);
L<theta> := ext<Qi|X^4 - (1 - 4*i)>;
u := (theta^3 - theta^2 - theta + 2*i - 1)/4;
s := hom<L->L|i*theta>;
pt := [1,u,s(u),s(s(u)),s(s(s(u))),u,s(u),s(s(u)),s(s(s(u)))];
ptcheck := C1(L)!pt;
assert Evaluate(D,pt) ne 0;
assert Evaluate(cc4,pt)^3 + 1728/2^7*Evaluate(D,pt)^10 eq 0;
assert Vector(pt) eq Vector(pt)*ChangeRing(M8,L)^4;

// A point on X(17) above j = 0

_<x> := PolynomialRing(K);
K3<z3> := ext<K|x^2+x+1>;
_<x> := PolynomialRing(K3);
L<c> := ext<K3|x^3 - sqrt17>;
ww := func<c,z3,z|(1/3*(-z^14 - z^12 + z^11 + 2*z^9 + 2*z^7 + z^6 + 3*z^5 +
  2*z^4 + 3*z^3 + 2*z^2 + 2*z + 1)*z3 + 1/3*(-2*z^15 - z^14 - 2*z^13 -
  z^12 + z^10 + 2*z^9 + z^8 + 2*z^7 + z^5 - z^4 + z^3 - z^2 + z -
  1))*c^2 + (1/3*(-z^15 - 2*z^14 - 2*z^12 + z^11 - z^10 + 2*z^9 +
  3*z^7 + z^6 + 4*z^5 + 2*z^4 + 4*z^3 + 3*z^2 + 2*z + 1)*z3 +
  1/3*(-3*z^15 - 2*z^14 - 3*z^13 - z^12 + z^10 + 3*z^9 + 2*z^8 + 3*z^7
  + 2*z^5 - 2*z^4 + z^3 - z^2 + z - 1))*c + 1/3*(-4*z^14 - 2*z^12 +
  2*z^11 - z^10 + 5*z^9 - z^8 + 5*z^7 + 2*z^6 + 6*z^5 + 4*z^4 + 8*z^3
  + 4*z^2 + 4*z + 2)*z3 + 1/3*(-4*z^15 - 3*z^14 - 4*z^13 - 2*z^12 +
  3*z^10 + 6*z^9 + 3*z^8 + 6*z^7 + 2*z^5 - 2*z^4 + 3*z^3 - 2*z^2 + 2*z
  - 2)>;
pt0 := [1] cat [ww(c,z3,z17^(9^i)) : i in [0..7]];
time check := C1(L)!pt0;
time assert Evaluate(cc4,pt0) eq 0;
time assert Evaluate(D,pt0) ne 0;

// Numerical approximations for the points above j = 0,1728,infty

CC<i> := ComplexField(100);
z17C := Exp(2*Pi(CC)*i/17);

err,nu := Minimum([Abs(z - z17C) : z in Conjugates(z17)]);
assert err lt 10^(-50);
M2C := Matrix(CC,9,9,[Conjugates(M2[i,j])[nu] : i,j in [1..9]]);
M8C := ChangeRing(M8,CC);
M17C := Matrix(CC,9,9,[Conjugates(M17[i,j])[nu] : i,j in [1..9]]);
apply := func<mat,vec|vec*Transpose(mat)>; // to give a left action
P<x> := PolynomialRing(QQ);
rep_pts := [];
all_pts := [];

for ctr in [1..3] do
  case ctr :
  when 1 : // j = infinity
    L := K;
    pt := [0,1,0,0,0,0,0,0,0];
  when 2 :  // j = 0
    L6<a> := NumberField(x^6 + x^4 - 11*x^2 + 12);
    L<z17> := ext<L6|CyclotomicPolynomial(17)>;
    Embed(K,L,z17);
    z3 := 1/4*(a^5 + 3*a^3 - 7*a - 2);
    c := -(2*a^4 + 7*a^2 - 13)/sqrt17;
    pt := [1] cat [ww(c,z3,z17^(9^i)) : i in [0..7]];
  when 3 : // j = 1728
    L<b> := NumberField(x^8 - 2*x^4 + 17);
    i := (1 - b^4)/4;
    s := hom<L->L|i*b>;
    u := (b^3 - b^2 - b + 2*i - 1)/4;
    pt := [1,u,s(u),s(s(u)),s(s(s(u))),u,s(u),s(s(u)),s(s(s(u)))]; 
  end case;
  check := C1(L)!pt;
  pt := Vector(L,pt);

  M0 := case<ctr| 1:M17, 2:M3, 3:M8^4, default:"">;
  zz := case<ctr| 1:z17, 2:z3, 3:1, default:"">;
  assert apply(ChangeRing(M0,L),pt) eq zz*pt;

  if ctr in [1,2] then
    err,nu := Minimum([Abs(z - z17C) : z in Conjugates(z17)]);
    assert err lt 10^(-50);
  else
    Phi17 := PolynomialRing(L)!CyclotomicPolynomial(17);
    assert IsIrreducible(Phi17);
    nu := 1;
  end if;
  pt := Vector([CC|Conjugates(x)[nu] : x in Eltseq(pt)]);
  
  case ctr :
  when 1 :
    gg := [M17^l : l in [0..16]];
    gg cat:= [M17^k*M2*g : k in [1..8],g in gg];
    gg := [M2^j*g : j in [0,1],g in gg];
    gg := [M8^i*g : i in [0..7],g in gg];
    assert #gg eq #G and #Set(gg) eq #G;
    pts := [pt] cat [apply(M17C^k*M2C,pt): k in [1..8]];
    pts := [apply(M2C^j,pt) : j in [0,1],pt in pts];
    pts := [apply(M8C^i,pt) : i in [0..7],pt in pts];
    assert #pts eq 144;
  when 2 : 
    gg := [M3^l : l in [0..2]];
    gg cat:= [M2*M17^k*g : k in [1,2,3,4,9],g in gg];
    gg := [M8^j*g : j in [0..7],g in gg];
    gg := [M17^i*g : i in [0..16],g in gg];
    assert #gg eq #G and #Set(gg) eq #G;
    pts := [pt] cat [apply(M2C*M17C^k,pt) : k in [1,2,3,4,9]];
    pts := [apply(M8C^j,pt) : j in [0..7],pt in pts];
    pts := [apply(M17C^i,pt) : i in [0..16],pt in pts];
    assert #pts eq 816;
  when 3 : 
    gg := [M8^k : k in [0..7]];
    gg cat:= [M2*M17^j*g : j in [0..16],g in gg];
    gg := [M17^i*g : i in [0..16],g in gg];
    assert #gg eq #G and #Set(gg) eq #G;
    pts := [apply(M8C^k,pt) : k in [0..3]];
    pts cat:= [apply(M2C*M17C^j,pt) : j in [0..16], pt in pts];
    pts := [apply(M17C^i,pt) : i in [0..16], pt in pts];
    assert #pts eq 1224;
  end case;
  
  rep_pts cat:= [Eltseq(pts[1])];
  all_pts cat:= [Eltseq(P): P in pts];
  
end for;

assert #rep_pts eq 3;
assert #all_pts eq 144 + 816 + 1224;

IndentPop();

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

printf "\nSection 5 : Computations in the symplectic case\n\n";
IndentPush();

print "Some random Fp-points on X(17)";

p := 10^9 + 7;
R9p := PolynomialRing(GF(p),9);
Cp := Curve(Proj(R9p),[R9p!f : f in X17eqns]);
ptsCp := [Cp!pt : pt in random_modp_points];
Rprint := RealField(8);

function IsBiInvariant(k,inv);
  assert k in [1,3];
  R9 := PolynomialRing(K,9);
  P8P8K := ProductProjectiveSpace(K,[8,8]);
  RRK := CoordinateRing(P8P8K);
  xx := [RRK.i : i in [1..9]];
  yy := [RRK.i : i in [10..18]];
  M2alt := k eq 1 select M2 else M2tilde;
  M17alt := k eq 1 select M17 else M17tilde;
  M2subst := [Evaluate(R9.i^M2,xx) : i in [1..9]];
  M2subst cat:= [Evaluate(R9.i^M2alt,yy) : i in [1..9]];
  M17subst := [Evaluate(R9.i^M17,xx) : i in [1..9]];
  M17subst cat:= [Evaluate(R9.i^M17alt,yy) : i in [1..9]];
  return (Evaluate(RRK!inv,M17subst) eq RRK!inv)
     and (Evaluate(RRK!inv,M2subst) eq RRK!inv);
end function;

for k in [1,3] do // the symplectic and antisymplectic cases

  if k eq 3 then
    IndentPop();
    printf "\nSection 6 : Computations in the anti-symplectic case\n\n";
    IndentPush();
  end if;

  P8P8 := ProductProjectiveSpace(QQ,[8,8]);
  RR<x0,x1,x2,x3,x4,x5,x6,x7,x8,
     y0,y1,y2,y3,y4,y5,y6,y7,y8> := CoordinateRing(P8P8);
  xx := [RR.i : i in [1..9]];
  yy := [RR.i : i in [10..18]];
  iota := k eq 1 select yy cat xx
                   else yy cat [-x0,-x2,-x3,-x4,-x5,-x6,-x7,-x8,-x1];

// Using the character table of G to compute the dimenision of the
// space of bi-invariants of degree (r,s)

  CT := CharacterTable(PSL(2,17));
  theta := CT[2];
  psi := k eq 1 select CT[2] else CT[3];
  dims := Matrix(6,6,[IP(Symm(theta,[r]),Symm(psi,[s])): r,s  in [0..5]]);
  assert dims eq (k eq 1 select Matrix(6,6,[
      [ 1,   0,   1,   1,   2,   2 ],
      [ 0,   1,   1,   2,   4,   7 ],
      [ 1,   1,   4,   6,  16,  30 ],
      [ 1,   2,   6,  16,  40,  94 ],
      [ 2,   4,  16,  40, 118, 277 ],
      [ 2,   7,  30,  94, 277, 694 ]])
    else Matrix(6,6,[
      [ 1,   0,   1,   1,   2,   2 ],
      [ 0,   0,   0,   1,   2,   6 ],
      [ 1,   0,   3,   5,  14,  29 ],
      [ 1,   1,   5,  15,  38,  93 ],
      [ 2,   2,  14,  38, 114, 275 ],
      [ 2,   6,  29,  93, 275, 693 ]]));

// Some bi-invariants obtained from the invariants Q,D and F.

  P<t> := PolynomialRing(RR);
  xy := [xx[i] + t*yy[i] : i in [1..9]];
  fij := func<f,d|[Coeff(Evaluate(f,xy),j) : j in [0..d]]>;
  Q20,Q11,Q02 := Explode(fij(Q,2));
  D30,D21,D12,D03 := Explode(fij(D,3));
  F40,F31,F22,F13,F04 := Explode(fij(F,4));
  
  print "The bi-invariants of degree (2,2)";

  if k eq 1 then
    HQinv := ChangeRing(HQinv,QQ);
    Dx := HQinv*H(D30);
    Dy := HQinv*H(D21);
    Fxy := HQinv*H(F31);
    AA := [ Q20*Q02, Q11^2, F22 ];
    AA cat:= [-2*AA[1] + AA[2] - (1/8)*Trace(Dx*Dy*Fxy)];
    AA0 := [Q11,AA[3],AA[4]];
  else
    M8M8 := DiagonalJoin(M8,M8);
    cyclic_sum := func<f|&+[f^(M8M8^j): j in [0..7]]>;
    AA := [ Q20*Q02, 4*x0^2*y0^2 + cyclic_sum(2*x0*x1*y3*y4
      + 2*x4*x5*y0*y1 + x1^2*y1*y7 + x2*x8*y1^2 + x1*x2*y2*y5
      + x3*x6*y1*y2 + x1*x3*y5*y8 + x1*x6*y1*y3 + (1/2)*x1*x5*y2*y6
      + (1/2)*x3*x7*y1*y5), 4*x0^2*y0^2 + cyclic_sum(2*x0*x1*y2*y8
      + 2*x1*x3*y0*y1 + x1^2*y2*y3 + x3*x4*y1^2 + x1*x2*y2*y5
      + x3*x6*y1*y2 + x1*x3*y5*y8 + x1*x6*y1*y3 + (1/2)*x1*x5*y1*y5
      + (1/2)*x2*x6*y1*y5) ];
    AA0 := [AA[2],AA[3]];
  end if;

  assert forall{f : f in AA | Evaluate(f,iota) eq f};
  assert forall{f : f in AA | IsBiInvariant(k,f) };
  mons := SetToSequence(&join[Set(Monomials(f)): f in AA]);
  mat := Matrix(QQ,#AA,#mons,[MC(f,mon): mon in mons,f in AA]);
  assert Rank(mat) eq #AA;

  if k eq 3 then
    print "The bi-invariants of degrees (3,1) and (1,3)";
    S31 := -16*x0^3*y0 + cyclic_sum((3*x0*x1*x5 + 3*x1^2*x4)*y0
      + (6*x0*x1*x3 + 6*x0*x4*x5 + 3*x2^2*x3 + 3*x1*x4^2 + x5^3
      + 3*x4*x6^2 + 6*x1*x2*x7 + 6*x3*x5*x8 + 6*x6*x7*x8)*y1);
    S13 := Evaluate(S31,iota);
    time assert IsBiInvariant(k,S31);
  end if;

  print "The bi-invariants of degree (3,3)";
  
  if k eq 1 then
    BB1 := [Q11*AA[2],Q11*AA[3],Q11*AA[4],D30*D03,D12*D21,Trace(Dx^3*Dy^3),
      Trace(Dx^2*Dy*Dx*Dy^2),Trace(Dx^2*Dy^2*Fxy),Trace(Dx*Dy*Dx*Dy*Fxy),
      Trace(Dx*Dy*Fxy^2),Trace(Fxy^3),Q11*AA[1],Q20*F13 + Q02*F31,
      &+[2*HQinv[i,j]*Deriv(F40,xx[i])*Deriv(F04,yy[j]): i,j in [1..9]] ];    
  else
    Hxx := func<f|HQinv*Matrix(9,9,
      [Deriv(Deriv(f,RR.i),RR.j): i,j in [1..9]])>;
    Hxy := func<f|HQinv*Matrix(9,9,
      [Deriv(Deriv(f,RR.i),RR.(j+9)): i,j in [1..9]])>;
    Hyx := func<f|HQinv*Matrix(9,9,
      [Deriv(Deriv(f,RR.(i+9)),RR.j): i,j in [1..9]])>;
    Hyy := func<f|HQinv*Matrix(9,9,
      [Deriv(Deriv(f,RR.(i+9)),RR.(j+9)): i,j in [1..9]])>;
    P := 2*(Q20*S13 + Q02*S31);
    Theta := func<i,j|Trace(Hxx(D30)*Hxy(AA[i])*Hyy(D03)*Hyx(AA[j]))>;
    Psi := func<i,j|Trace(Hxx(S31)*Hxy(AA[i])*Hyx(AA[j]))>;
    U2 := Trace(Hxy(P)*Hyx(AA[2]));
    U3 := Trace(Hxy(P)*Hyx(AA[3]));
    V := Trace(Hxy(Theta(1,2))*Hyx(AA[3]));
    BB1 := [D30*D03,P,Theta(1,2),Theta(1,3),Theta(2,2),
      Theta(2,3),Theta(3,3),U2,U3,Psi(2,2) + Evaluate(Psi(2,2),iota),
      Psi(2,3) + Evaluate(Psi(2,3),iota),V];
  end if;

  assert forall{B : B in BB1 | Evaluate(B,iota) eq B};

// The change of basis matrices recorded in Appendix A

  cob := k eq 1 select Matrix(QQ,14,14,[
    [ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ],
    [ 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ],
    [ 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ],
    [ -24, -12, 28, 8, -4, 0, 4, 4, 0, -4, 4, 0, 0, 0 ],
    [ -34, -15, 34, 7, -8, -3, 0, 5, -2, -9, 7, 0, -1, 0 ],
    [ -480, -288, 320, 160, -64, -32, 192, 160, 0, -96, 96, 512, -96, 0 ],
    [ 248, 144, -96, -144, 0, -16, -192, -144, -32, 48, 16, 448, -48, 16 ],
    [ -24, -56, -128, -40, 0, -24, 32, 56, -16, -24, -8, -32, -40, 8 ],
    [ -384, -144, 688, 96, -64, 0, 32, -48, -16, 0, 128, -64, 0, 8 ],
    [ 328, 204, -240, -76, 32, 4, -144, -92, 8, 36, -20, 24, 20, -12 ],
    [ -366, -180, 96, 72, -48, 12, 192, 168, -12, -48, -84, 36, 0, 36 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 ] ])
  else Matrix(QQ,12,12,[
    [ -4, -12, 8, -4, -8, -8, 0, -4, 0, 0, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0 ],
    [ -56, -176, 148, -68, -160, -140, 36, -52, 24, 24, -8, 4 ],
    [ -104, -240, 156, -76, -128, -164, 12, -92, 8, 8, 8, 4 ],
    [ -16, -128, 124, -36, -144, -92, 12, -36, 8, 8, 8, 12 ],
    [ -40, -48, 48, -24, -32, -40, 48, -24, 32, 32, 0, 16 ],
    [ -64, -176, 132, -60, -128, -148, 4, -76, -8, -8, 24, -12 ],
    [ -36, -72, 96, -96, 12, 24, -36, -12, 0, 24, 24, -12 ],
    [ -60, -48, 48, -24, -12, -48, 60, -60, 48, 48, 72, 24 ],
    [ 36, 72, -72, 24, 60, 24, -60, 12, -72, -72, 24, -48 ],
    [ -72, -84, 108, -48, 0, 24, -12, -24, 12, 24, 0, 12 ],
    [ -208, 204, -72, 28, 364, 116, 288, -44, 208, 224, 32, 104 ] ]);

  cobinv := cob^(-1);
  BB := [&+[cobinv[i,j]*BB1[j] : j in [1..#BB1]]: i in [1..#BB1]];
  assert BB1 eq [&+[cob[i,j]*BB[j] : j in [1..#BB]]: i in [1..#BB]];

  printf "Creating the surface Sigma_%o\n",k;
  
  if k eq 1 then 
    R<z1,z2,z3,z4,z5,z6,z7,z8,z9,z10,z11> := PolynomialRing(QQ,11);
  else 
    R<z1,z2,z3,z4,z5,z6,z7,z8,z9,z10> := PolynomialRing(QQ,10);
  end if;
  quads,cubics := SigmakEquations(k,R);
  assert #quads eq (k eq 1 select 13 else 14);
  assert #cubics eq (k eq 1 select 21 else 2);
  Sigmak := Scheme(Proj(R),quads cat cubics);
  assert Dimension(Sigmak) eq 2;
  assert Degree(Sigmak) eq (k eq 1 select 29 else 24);
  time assert {* Degree(f): f in MinimalBasis(Sigmak) *}
    eq {* 2^^#quads, 3^^#cubics *};

  // reln8 := u*z5*z10 - v*(z1*z5 - z3*z5 + z1*z6 - z5*z8 + z5*z9);
  // reln9 := u*v^2*z7 - (z2 - z3 - z5)*(z6 - z7)*z10 + z5*z8*z10;
  if k eq 3 then
    p1 := z5*z10;
    p2 := z1*z5 - z3*z5 + z1*z6 - z5*z8 + z5*z9;
    q1 := z7;
    q2 := (z2 - z3 - z5)*(z6 - z7)*z10 - z5*z8*z10;
  end if;

  print "Checking the relations";
  
  for method in ["modp","numerical"] do
    print "Method =",method;
    iszero := method eq "modp" select func<x|x eq 0>
                                 else func<x|Abs(x) lt 10^(-80)>;
    for ctr in [1..100] do // more iterations => more check
      if method eq "modp" then
        pt1 := Eltseq(Random(ptsCp));
        pt2 := Eltseq(Random(ptsCp));
      else
        pt1 := Random(rep_pts);
        pt2 := Random(all_pts);
      end if;
      bb := [Evaluate(B,pt1 cat pt2): B in BB];
      assert forall{ i : i in [Rank(R)+1..#BB] | iszero(bb[i]) };
      xx := [bb[i] : i in [1..Rank(R)]];
      assert forall{r : r in quads cat cubics | iszero(Evaluate(r,xx))};
      if k eq 3 then
        aa := [Evaluate(AA[i],pt1 cat pt2): i in [2..3]];
        diff1 := aa[1]*Evaluate(p1,xx) - aa[2]*Evaluate(p2,xx);
        diff2 := aa[1]*aa[2]^2*Evaluate(q1,xx) - Evaluate(q2,xx);
        assert iszero(diff1) and iszero(diff2);
      end if;	
    end for;
  end for;

  if k eq 1 then 
    print "Proposition 5.3 : Showing Sigma_1 is birational to Sigma";
  else
    print "Proposition 6.2 : Showing Sigma_3 is birational to Sigma";
  end if;

  qq := QuadricIntersection(k,z5/z6,z1,z2,z3,z4);
  assert [Numerator(q): q in qq] subset Ideal(Sigmak);
  R4<z1,z2,z3,z4> := PolynomialRing(QQT,4);
  C4 := Curve(Proj(R4),QuadricIntersection(k,T,z1,z2,z3,z4));
  if k eq 1 then
    uu := [ T*(T + 1)*z1, T*(z1 + z2 - 2*z3 - T*z4),
            (T + 1)*z3, T*(-z1 + 2*z3 + T*z4) ];
  else
    uu := [ z1 - T*z2, (T + 2)*z1 + z2 - (T + 1)*(z3 - z4),
            (T + 1)*z3, (T + 1)*z4 ];
  end if;
  R4<u1,u2,u3,u4> := R4;
  C4a := Curve(Proj(R4),QI_Simplified(k,T,u1,u2,u3,u4));
  gamma := map<C4->C4a|uu>;
  if k eq 1 then 
    psi := map<C4a->E|[-T*u1/u4,T*u1*(u1 - u2 - u4)/(u2*u4),1]>;
  else
    psi := map<C4a->E|[T*(u1 + T*u2)/u4,T*(u1 + T*u2)^2/(u1*u4),1]>;
  end if;
  KE<x,y> := FunctionField(E);
  if k eq 1 then 
    psiinv := map<E->C4a|[ x*(T*y - (2*T + 1)*x + T^2*(T + 1)),
      -T*x^2 + (T + 1)*y + (T^3 - T^2 - 2*T - 1)*x + T^3*(T + 1),
      T*x^2 - (T + 1)^2*y - (T - 2)*(T + 1)^3*x - T^3*(T + 1)^2,
      -T*(T*y - (2*T + 1)*x + T^2*(T + 1))]>;
  else
    psiinv := map<E->C4a|[ T*x^2*(x + T), x*(x + T)*(y - x),
      -T*x*(y + x^2 - x), T^2*(x + T)*y ]>;
  end if;
  assert psiinv*psi eq IdentityMap(E);
  assert psi*psiinv eq IdentityMap(C4a);

  print "Computing the inverse";

  // Solving for z1,z2,z3,z4
  
  zz := DefiningEquations(psiinv*gamma^(-1));
  zz := [Evaluate(z,[x,y,1]): z in zz];
  
  if k eq 1 then
    zz := [T^2*(T + 1)*f : f in zz];
    assert zz eq [
      T*(T*x*y - (2*T + 1)*x^2 + T^2*(T + 1)*x),
      -T*(T + 1)*(T*x^2 + (T^2 - T - 1)*y - (T - 1)*(T + 1)^2*x),
      T^3*x^2 - T^2*(T + 1)^2*(y + (T + 1)*(T - 2)*x + T^3),
      T*x*y - (2*T^2 + 2*T + 1)*x^2 + T*(T + 1)*(T + 2)*y 
        + T*(T + 1)*(2*T^3 - 3*T - 3)*x + T^3*(T + 1)^2*(2*T - 1) ];
  else
    zz := [(T + 1)^2*z : z in zz];
    assert zz eq [
      T*(x^2 - T^2*x - T^3)*y - T^2*(x^3 - x^2),
      (x^2 - T^2*x - T^3)*y - (T^2 + 3*T + 1)*x^3 - T^2*(T + 2)*x^2,
      -T*(T + 1)*(x*y + x^3 - x^2),
      T^2*(T + 1)*(x + T)*y ];
  end if;

  // Some polynomials to help solve for the remaining z_i
  
  if k eq 1 then
    R<z1,z2,z3,z4,z5,z6,z7,z8,z9,z10,z11> := R;
    function nd(T,z1,z2,z3,z4)
      num := (T + 4)*z1^2 + T*(T + 4)*z1*z2 + 2*(T^2 + 2*T - 2)*z1*z3
        + T*(T^2 - 2*T + 3)*z1*z4 + T*(T + 1)*z2^2 - 3*T*z2*z3 - T^2*z2*z4
        - 2*T*(T + 2)*z3^2 - 3*T*(T^2 + 2)*z3*z4 - T^2*(T^2 - T + 3)*z4^2;
      den := (T - 2)*z1 - 2*(T - 1)*z3 - T*(T - 1)*z4;
      return num,den;
    end function;
    nums := [ (-3*z1^2 - z1*z2 + 5*z1*z3 + z1*z4 - z1*z5 - z2*z3)*z5,
      (3*z1^2 - 6*z1*z3 + z1*z5 + z1*z7 - z2^2 + z2*z3)*z6, -7*z1^2*z5
      - 11*z1*z2*z5 - z1*z2*z6 + 7*z1*z3*z5 - 4*z1*z3*z6 + z1*z4*z5 -
      z1*z4*z6 - 2*z1*z5^2 + z1*z5*z6 - z1*z5*z7 + 7*z1*z5*z8 +
      2*z1*z6*z8 - 4*z2^2*z5 + 11*z2*z3*z5 + z2*z3*z6 + 2*z2*z4*z5 -
      2*z2*z5^2 - 2*z2*z5*z7 + 4*z2*z5*z8 - 2*z2*z6*z7 + 10*z3^2*z5 +
      4*z3^2*z6 + z3*z4*z5 + z3*z4*z6 - 2*z3*z5^2 - z3*z5*z6 -
      8*z3*z5*z7 - 11*z3*z5*z8 - 2*z3*z6*z7 - 2*z3*z6*z8 - 2*z4*z5*z8
      + z5^2*z7 + 2*z5^2*z8 + 2*z5*z7^2 + 3*z5*z7*z8 + 2*z6*z7*z8,
      -14*z1^2 + 32*z1*z2 - 200*z1*z3 - 176*z1*z4 + 18*z1*z5 -
      51*z1*z6 + 22*z1*z7 + 21*z1*z8 - 12*z1*z9 + 42*z2^2 - 112*z2*z3
      - 39*z2*z4 + 20*z2*z5 - 5*z2*z6 - 7*z2*z7 - 28*z2*z8 + 120*z3^2
      + 288*z3*z4 - 24*z3*z5 + 59*z3*z6 + 60*z3*z7 + 55*z3*z8 +
      21*z3*z9 + 57*z4^2 - 69*z4*z5 + 19*z4*z6 - 35*z4*z7 - 18*z4*z8 -
      18*z5*z6 - 15*z5*z7 - 6*z5*z8 - 3*z5*z9 - 14*z6*z7 - 2*z6*z8 -
      2*z7^2 + 7*z7*z8, 40*z1^2 + 13*z1*z2 - 171*z1*z3 - 19*z1*z4 +
      12*z1*z5 - 6*z1*z6 + 15*z1*z7 - 3*z1*z8 - 11*z1*z9 - z2^2 -
      33*z2*z3 + 4*z2*z4 - z2*z5 + z2*z7 + 183*z3^2 + 50*z3*z4 -
      30*z3*z5 + 12*z3*z6 - 40*z3*z7 + z3*z8 + 19*z3*z9 + 3*z4^2 -
      5*z4*z5 + z4*z6 - 10*z4*z7 - 7*z4*z8 + 3*z4*z9 + 3*z5*z7 + z5*z8
      - 3*z5*z9 + z6*z9 + 6*z7^2 + 5*z7*z8 + z8^2 - z8*z9 ];
    dens := [ z1*(z5 + z6), z1*z6 - 2*z3*z6 - z4*z5, (2*z1 - 5*z3 + z5
      + 2*z7)*z5, 15*z1 - 14*z2 + 8*z3 + 24*z4 + 8*z6 + 8*z7, 31*z1 +
      9*z2 - 66*z3 - 15*z4 + 12*z5 - 5*z6 + 13*z7 + 4*z8 ];
  else
    R<z1,z2,z3,z4,z5,z6,z7,z8,z9,z10> := R;
    nums := [ z1*z4 + z2^2 - z2*z3 + 2*z2*z4 + z2*z5 - z3*z4 - z3*z5 +
      z4^2 + z4*z5, z1*z6 - z2*z5, z1*z2 + z1*z7 + 2*z2^2 - z2*z3 +
      z2*z4 + z2*z5 + 2*z2*z6 + z4*z6 + z5*z6 + z6^2 - z6*z7, z1*z2 +
      z3*z5 - z3*z7 + z4*z6 - z4*z7 - z5*z8, z3*z7 - z4*z6 + z4*z7 +
      z5*z8 + z6*z9 ];
    dens := [ z1, -z1 - z2 - z4, -z2 - z4, z2 - z3 + z4 + z6, -z5 ];
  end if;

  // Checking these polynomials vanish on Sigmak

  if k eq 1 then 
    num,den := nd(z5/z6,z1,z2,z3,z4);
    assert Numerator(num - z5*den) in Ideal(Sigmak);
  end if;
  m := Rank(R);  
  for i in [1..5] do
    assert forall{k : k in [m-5+i..m] | Degree(nums[i],k) eq 0};
    assert forall{k : k in [m-5+i..m] | Degree(dens[i],k) eq 0};
    assert nums[i] - R.(m-5+i)*dens[i] in Ideal(Sigmak);
  end for;

  // Solving for the remaining z_i
  
  if k eq 1 then
    zz := [z/zz[1] : z in zz];
    num,den := nd(T,zz[1],zz[2],zz[3],zz[4]);
    zz cat:= [num/den,num/den/T,0,0,0,0,0];
  else 
    P<X> := PolynomialRing(KE);
    zzz := [P!f : f in zz] cat [T*X,X,0,0,0,0];
    poly := Evaluate(nums[1] - z6*dens[1],zzz);
    rt := Roots(poly)[1,1];
    zz cat:= [T*rt,rt,0,0,0,0];
    nums := [nums[i] : i in [2..5]];
    dens := [dens[i] : i in [2..5]];
  end if;
  
  for i in [1..m-6] do
    zz[i+6] := Evaluate(nums[i]/dens[i],zz);
  end for;

  if k eq 3 then
    la := Evaluate(p2/p1,zz);
    mu := Evaluate(q2/q1,zz);
    zz := [la^2*mu,la*mu] cat [la^2*mu*f : f in zz];
    wts := [2,2,3,3,3,3,3,3,3,3,3,3];
    chi := (x - T^2)*(T*y - x^2 - T*(T+2)*x)/(T^2*(T + 1));
    zz := [zz[i]/chi^wts[i] : i in [1..12]];
  end if;

  // Checking against the answer recorded in "congr17-formulae.m"

  assert zz eq InverseMap(k,T,x,y);

  // Checking our formula for the inverse does indeed have
  // image contained in Sigmak
 
  zz1 := (k eq 1 select zz else [zz[i+2] : i in [1..10]]);
  time assert forall{f : f in quads | Evaluate(f,zz1) eq 0};
  time assert forall{f : f in cubics | Evaluate(f,zz1) eq 0};
  if k eq 3 then 
    assert zz[1]*Evaluate(p1,zz1) eq zz[2]*Evaluate(p2,zz1);
    assert zz[1]*zz[2]^2*Evaluate(q1,zz1) eq Evaluate(q2,zz1);
  end if;

  print "Extra relations for computing the j-maps";

  if k eq 1 then
    S := PolynomialRing(QQ,[1,2,2,3,3,3,3,3,3,3,3]);
    gg,hh := ExtraRelationsforJmap1(S);
    assert [Degree(g): g in gg] eq [ 3, 3, 9, 15, 15 ];
    assert [Degree(h): h in hh] eq [ 0, 0, 3, 5, 4 ];
  else
    S := PolynomialRing(QQ,[2,2,3,3,3,3,3,3,3,3,3,3]);
    gg,hh := ExtraRelationsforJmap3(S);
    assert [Degree(g): g in gg] eq [ 3, 6, 8, 8, 12, 18, 18 ];
    assert [Degree(h): h in hh] eq [ 0, 2, 3, 2, 5, 8, 6 ];
  end if;

  // Matrices to relate to an old choice of basis.
  // (In principle these could be made redundant by updating
  // gg and hh to use the new bases - however some care
  // would be needed to avoid slowing down the verificaton
  // that follows.)
  N := k eq 1 select Matrix(QQ,8,11,[
    [ 24, 12, -28, -8, 4, 0, -4, -4, 0, 4, -4 ],
    [ -4, -1, 5, 2, -1, 0, 0, 0, 0, -1, 1 ],
    [ -4, -2, 4, 1, -1, 0, 0, 1, 0, -1, 0 ],
    [ 13, 6, -16, -5, 3, 0, 0, -2, 0, 3, -2 ],
    [ -6, -3, 4, 2, -1, 0, 2, 2, 0, -1, 0 ],
    [ 8, 2, -12, -4, 2, -1, 0, 0, 0, 1, 0 ],
    [ -5, -4, 1, 2, -1, 0, 2, 3, 0, -2, -1 ],
    [ -15, -7, 19, 3, -4, -1, -2, 1, -1, -4, 3 ]])
  else Matrix(QQ,10,10,[
    [ 0, -1, 1, 0, -1, -1, 0, 0, 0, 0 ],
    [ 1, 2, -1, 1, 1, 1, 0, 0, 0, 0 ],
    [ -1, -1, 1, -1, 0, 0, 0, 0, 0, 0 ],
    [ 1, 1, -1, 0, 0, 0, 0, 0, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, -1, 1, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 ],
    [ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 ],
    [ 0, 0, 0, 0, -1, -1, 0, 0, 0, 0 ]]);

  if k eq 3 then
    print "The bi-invariants of degrees (3,2) and (2,3)";
    Chi := func<i,j|Trace(Hxx(D30)*Hxy(AA[i])*Hyx(AA[j]))>;
    T32 := -1/12*(4*Chi(1,1) - 3*Chi(1,3) + 3*Chi(2,3));
    T23 := Evaluate(T32,iota);
  end if;

  for method in ["modp","numerical"] do
    print "Method =",method;
    iszero := method eq "modp" select func<x|x eq 0>
                                 else func<x|Abs(x) lt 10^(-50)>;
    for ctr in [1..100] do // more iterations => more check
      if method eq "modp" then
        pt1 := Eltseq(Random(ptsCp));
        pt2 := Eltseq(Random(ptsCp));
      else
        pt1 := Random(rep_pts);
        pt2 := Random(all_pts);
      end if;
      pt := pt1 cat pt2;
      aa := [Evaluate(f,pt): f in AA0];
      bb := [Evaluate(f,pt): f in BB];
      xx := aa cat [&+[N[i,j]*bb[j]: j in [1..Ncols(N)]]: i in [1..Nrows(N)]];
      if k eq 1 then
        alpha := [ Evaluate(D30,pt)*Evaluate(D03,pt),
          Evaluate(D21,pt)*Evaluate(D12,pt),
          Evaluate(D21,pt)^3*Evaluate(D03,pt) 
            + Evaluate(D12,pt)^3*Evaluate(D30,pt), 
          Evaluate(cc4,pt1)*Evaluate(cc4,pt2),
          Evaluate(cc4,pt1)*Evaluate(D03,pt)^3*Evaluate(D12,pt)
            + Evaluate(cc4,pt2)*Evaluate(D30,pt)^3*Evaluate(D21,pt) ];
      else
        alpha := [ Evaluate(D30,pt)*Evaluate(D03,pt),
          Evaluate(S31,pt)*Evaluate(S13,pt),
          Evaluate(T32,pt)*Evaluate(T23,pt),
          Evaluate(D03,pt)*Evaluate(S31,pt)*Evaluate(T32,pt)
            + Evaluate(D30,pt)*Evaluate(S13,pt)*Evaluate(T23,pt),
          Evaluate(S31,pt)*Evaluate(T23,pt)^2
            + Evaluate(S13,pt)*Evaluate(T32,pt)^2,	
          Evaluate(cc4,pt1)*Evaluate(cc4,pt2),
          Evaluate(D03,pt)^3*Evaluate(T23,pt)*Evaluate(cc4,pt1)
            + Evaluate(D30,pt)^3*Evaluate(T32,pt)*Evaluate(cc4,pt2) ];
      end if;
      assert forall{i : i in [1..#alpha] |
        iszero(Evaluate(gg[i],xx) - alpha[i]*Evaluate(hh[i],xx))};
    end for;
  end for;

  if k eq 1 then
    zz := [1,zz[2],zz[3]] cat [&+[N[i,j]*zz[j]: j in [1..11]]: i in [1..8]];
  else
    mat := DiagonalJoin(IdentityMatrix(QQ,2),N);
    zz := [&+[mat[i,j]*zz[j] : j in [1..12]] : i in [1..12]];
  end if;

  print "Solving for j1 + j2 and j1*j2";

  time alpha := [Evaluate(gg[i]/hh[i],zz) : i in [1..#gg]];

  print "Alternative way to represent the function field Q(Sigma)";
  K<x> := FunctionField(QQT);
  _<Y> := PolynomialRing(K);
  KE<y> := ext<K|Y^2 + (T + 1)*(T - 2)*x*Y + T^3*Y - (x^3 - x^2)>;
  assert Evaluate(Equation(E),[x,y,1]) eq 0;
  aa := eval Sprint(alpha);

  P<X> := PolynomialRing(KE);
  ff := k eq 1 select [
    X^2 - aa[3]*X + aa[1]*aa[2]^3,
    X^2 - aa[5]*X + aa[1]^3*aa[2]*aa[4]]
  else [
    X^2 - aa[4]*X + aa[1]*aa[2]*aa[3],
    X^2 - aa[5]*X + aa[2]*aa[3]^2,
    X^2 - aa[7]*X + aa[1]^3*aa[3]*aa[6]];
  FF := Z17Equation(k,T,x,y);
  mm := MySquareRoots(k,T,x,y);
  assert forall{ i : i in [1..#ff] | Discriminant(ff[i]) eq mm[i]^2*FF };
  L<c> := quo<P|X^2 - FF>;
  rr := [(-Coeff(ff[i],1) + mm[i]*c)/2 : i in [1..#ff]];
  if k eq 1 then 
    time jsum := -2^7*Trace(rr[2]^3/rr[1])/aa[1]^9; 
    time jprod := 2^14*aa[4]^3/aa[1]^10; 
  else
    time jsum := -2^7*aa[2]*Trace(rr[3]^3/(rr[1]*rr[2]))/aa[1]^9;
    time jprod := 2^14*aa[6]^3/aa[1]^10;
  end if;

  print "Checking j1 + j2 and j1*j2 against the stored answers";
  
  AA,BB,D1,D17 := Z17ModuliInformation(k,T,x,y);
  DD := D1*D17^17;
  time assert jprod eq AA^3/DD;
  time assert jprod - 1728*jsum + 1728^2 eq BB^2/DD;

  if k eq 1 then // Remark 5.4
    R4<z1,z2,z3,z4> := PolynomialRing(QQT,4);
    qq := QuadricIntersection(k,T,z1,z2,z3,z4);
    eqn1 := (1/T^2)*Resultant(qq[1],qq[2],z4);
    assert eqn1 eq
      T*(T + 1)^3*z1^3*z2 + T*(T + 1)*(T^2 + T + 1)*z1^2*z2^2
      + (T + 1)^4*z1^3*z3 + T*(T + 1)^2*(T^2 - T - 3)*z1^2*z2*z3
      + T^2*(T + 1)*(T^2 - 2*T - 1)*z1*z2^2*z3 - T^4*z2^3*z3
      + (T + 1)^4*(T - 1)*z1^2*z3^2 - T^2*(T + 1)^2*z1*z2*z3^2;
  end if;
  
end for;

IndentPop();

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

printf "\nSection 7 : Some modular curves\n\n";

QQt<t> := FunctionField(QQ);
PP<e> := LaurentSeriesRing(QQt,50);
R<X> := PolynomialRing(PP);

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

function BetterApproximation(pt)
  T,u,v := Explode(pt);
  cc := v^2 + (T^2 - T - 2)*u*v + T^3*v - u^3 + u^2;
  rts := [r : r in Roots(cc) | Valuation(r[1]) ge 0];
  assert #rts eq 1;
  return [Evaluate(f,rts[1,1]): f in pt];
end function;

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

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

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

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

for ctr in [1..#mm] do
  m,k,pt,c,r := Explode(Table1[ctr]);
  assert m eq mm[ctr];
  g := GenusX0N(m);
  assert g in [0,1];
  assert k eq (LegendreSymbol(m,17) eq 1 select 1 else 3);
  prec := 20;
  if m in [3,6] then prec := 40; end if; 
  if m eq 9 then prec := 30; end if;
  PP<e> := LaurentSeriesRing(QQt,prec);
  R<X> := PolynomialRing(PP);
  pt := eval Sprint(pt);
  
  printf "%o ",m;
  
  if r eq 0 then
    T,x,y := Explode(pt);
    assert Z17Equation(k,T,x,y) eq c;
  else
    T,x,y := Explode(BetterApproximation(pt));
    assert Z17Equation(k,T,x,y) + O(e^(r+1)) eq c*e^r + O(e^(r+1));
  end if;

  if r eq 0 then
    E := EC([ t^2 - t - 2, -1, t^3, 0, 0 ]);
    x1 := Coeff(x,0);
    y1 := Coeff(y,0);
    P := E![x1,y1];
    KE<x,y> := FunctionField(E);
    AA,BB,D1,D17 := Z17ModuliInformation(k,t,x,y);
    AA := Evaluate(AA,P);
    BB := Evaluate(BB,P);
    D1 := Evaluate(D1,P);
    D17 := Evaluate(D17,P);
  else
    AA,BB,D1,D17 := Z17ModuliInformation(k,T,x,y);
  end if;

  DD := D1*D17^17;
  jprod := AA^3/DD;
  jsum := (jprod + 1728^2 - BB^2/DD)/1728;

  if r ne 0 then
    jsum := Evaluate(jsum,0);
    jprod := Evaluate(jprod,0);
  end if;

  X0m := SmallModularCurve(m);
  
  if g eq 0 then
    KK := QQt;
    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(jsum,t1)*X + Evaluate(jprod,t1);
  assert {r[1] : r in Roots(jpoly)} eq {j,j1};

  if m eq 18 then // example in the text
    assert j1 eq Evaluate(j,6/t);
    assert t1 eq -t/((t+2)*(t+3));
    num := (t + 2)^12 - 8*(t + 2)^9 + 16*(t + 2)^3 + 16;
    den := (t + 2)^9*((t + 2)^3 - 8)*((t + 2)^3 + 1)^2;    
    assert j eq num^3/den;
  end if;

end for;

print "";

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

printf "\nSection 8 : Examples and further questions\n\n";
IndentPush();

Table2 := [
  < 1, [ 1/3, -2/75, -11/125 ], "279809270*", "3077901970*", 0 >,
  < 1, [ -3, -27, 108 ],        "1849a1", "1849a2", 43 >,
  < 1, [ -5/6, -5/24, 5/16 ],   "4489a1", "4489a2", 67 >,
  < 3, [ 1, 2, 4 ],             "27a2", "27a4", 27 >,
  < 3, [ 9/7, 27/49, -54/49 ],  "3675b1", "47775b1", 0 >,
  < 3, [ -5/14, 125/392, 375/1568 ], "1225h1", "1225h2", 37 >,
  < 3, [ 11/39, 1771/6591, 116380/257049 ], "26569a1", "26569a2", 163 >
];

EllCrv := func<ref|case<ref|
  "279809270*":E1, "3077901970*":E2, default:EC(ref) >>;

for dd in Table2 do

  k,pt,ref1,ref2,deg := Explode(dd);
  EE1 := EllCrv(ref1);
  EE2 := EllCrv(ref2);
  if deg ne 0 then
    flag,isog := IsIsogenous(EE1,EE2);
    assert flag;
    assert Degree(isog) eq deg;
  else
    assert not IsIsogenous(EE1,EE2);
  end if;
  aa := TracesOfFrobenius;
  vec1 := Vector(GF(17),aa(EE1,10^4));
  vec2 := Vector(GF(17),aa(EE2,10^4));
  assert #[x : x in Eltseq(vec1 - vec2) | x ne 0] le 1;
  
  T,x,y := Explode(pt);
  AA,BB,D1,D17 := Z17ModuliInformation(k,T,x,y);
  DD := D1*D17^17;
  jprod := AA^3/DD;
  jsum := (jprod + 1728^2 - BB^2/DD)/1728;
  Qt<t> := PolynomialRing(QQ);
  jpoly := t^2 - jsum*t + jprod;
  jj := {r[1] : r in Roots(jpoly)};
  assert jj eq {jInvariant(EE1),jInvariant(EE2)};
  
end for;

QQT<T> := QQT;
KE<x,y> := FunctionField(E);

Table_j_inf := [
< 1, [ 0, -T^3 ], T^16*(T - 1)^2*(T + 1)^4, [0,0,5,0] >,
< 1, [ T^3 + T^2, T^4 + 2*T^3 + T^2 ],
     T^18*(T + 1)^4*(T + 2)^4, [0,0,3,0] >,
< 1, [ T^3 - T, T^4 - T^2 - T ],
     -T^16*(T + 1)^4*(T - 1)*(3*T + 1)^3, [0,0,2,0] >,
< 3, [ -T, -T ], T^12*(T + 2)^4*(T^2 + T + 1)^4, [0,0,6,0] >,
< 3, [ -T, -T^2 - T ], T^12*(2*T + 1)^4*(T + 1)^8, [6,9,5,1] >,
< 3, [ T^2, -T^4 + T^2 ], T^18*(T + 1)^10, [0,0,5,0] >];

Table_spurious := [
< 1, [ -T, -T ], 4, T^8*(T^5 - 5*T^3 - 20*T^2 - 24*T - 19)
                   *(T^5 + 4*T^4 + 7*T^3 + 8*T^2 + 4*T + 1) >,
< 1, [ T^2, -T^4 + T^2 ], 2, -T^8*(T^3 + 2*T + 1)
     *(11*T^9 + 24*T^8 + 46*T^7 + 61*T^6 + 60*T^5 + 44*T^4
     + 21*T^3 + 4*T^2 - 2*T - 1) >,
< 3, [ T^3 + T^2, -T^5 - T^4 + T^2], 4, T^6*hh/(T^2 + 2*T + 2)^4 > ]
  where hh is 49*T^20 + 854*T^19 + 6959*T^18 + 35510*T^17 +
    127907*T^16 + 346802*T^15 + 735814*T^14 + 1252442*T^13 +
    1738217*T^12 + 1987496*T^11 + 1883412*T^10 + 1482532*T^9 +
    968144*T^8 + 521836*T^7 + 229843*T^6 + 81334*T^5 + 22489*T^4 +
    4638*T^3 + 654*T^2 + 52*T + 1;

for dd in Table_j_inf do
  k,P,h,v := Explode(dd);
  Fk := Z17Equation(k,T,x,y);
  AA,BB,D1,D17 := Z17ModuliInformation(k,T,x,y);
  assert Evaluate(Fk,E!P) eq h;
  assert [Valuation(f,E!P) : f in [AA,BB,D1,D17]] eq v;
end for;

for dd in Table_spurious do
  k,P,v,h := Explode(dd);
  Fk := Z17Equation(k,T,x,y);
  assert Evaluate(Fk/(x - P[1])^v,E!P) eq h;
end for;

// The curves listed in the statement of Conjecture 8.2

curvelist := [[
    [ 0, -T^3 ],                        // at infinity
    [ -T, -T^2 - T ],                   // corresponds to X0(18)
    [ -T, -T],                          // spurious
    [ T^2, T^2 ],                       // corresponds to X0(25)
    [ T^2, -T^4 + T^2 ],                // spurious
    [ T^3 + T^2, T^4 + 2*T^3 + T^2 ],   // at infinity
    [ T^3 - T, T^4 - T^2 - T]],         // at infinity
    [
    [ -T, -T ],                         // at infinity
    [ -T, -T^2 - T ],                   // at infinity
    [ T^2, -T^4 + T^2 ],                // at infinity
    [ T^3 + T^2, -T^5 - T^4 + T^2]]];   // spurious

for k in [1,3] do
  curves := k eq 1 select curvelist[1] else curvelist[2];
  printf "Some curves on Z(17,%o)\n",k;
  for P in curves do
    mm := Sprint([Sprint(a): a in P]);
    printf "(x,y) = %o \t",mm cat &cat[" ": i in [#mm..32]];
    if exists{d : d in Table_j_inf | d[1] eq k and d[2] eq P} then
      print " at infinity";
    end if;
    if exists{d : d in Table_spurious | d[1] eq k and d[2] eq P} then
      print " spurious";
    end if;
    if exists(d){d : d in Table1 | d[2] eq k and d[3] eq [T] cat P} then
      printf " corresponds to X0(%o)\n",d[1];
    end if;
  end for;
end for;

print "Searching for rational points on Z(17,1) and Z(17,3)";

P<u> := PolynomialRing(QQT);
PP<v> := PolynomialRing(P);
Weqn := Evaluate(Equation(E),[u,v,1]);
cc := [Coefficient(Discriminant(Weqn),i): i in [0..3]];
PQ := PolynomialRing(Rationals());
PZ := PolynomialRing(Integers());

for b in [1..50] do
  for a in [1..b] do
    if GCD(a,b) eq 1 then
      for T in {a/b,b/a,-a/b,-b/a} do
        g := PQ![Evaluate(c,T): c in cc];
        den := LCM([Denominator(c): c in Coefficients(g)]);
	g := PZ!(den^2*g);
	gg,c0 := Factorization(g);
	c0 := Integers()!PowerFreePart(Rationals()!c0,2);
	h := c0*(&*[g[1] : g in gg | IsOdd(g[2])]);
	assert IsSquare(g/h);
   	C := HyperellipticCurve(h);	
        pts := Points(C : Bound := 10^4);
  	xlist := {P[1]/P[3] : P in pts | P[3] ne 0};
  	for x in xlist do
  	  f1 := PQ.1^2 + ((T^2 - T - 2)*x + T^3)*PQ.1 - x^3 + x^2;
	  ylist := [r[1] : r in Roots(f1)];
	  for y in ylist do
	    assert y^2 + (T^2 - T - 2)*x*y + T^3*y eq x^3 - x^2;
	    for k in [1,3] do
              Fk := Z17Equation(k,T,x,y);
	      curves := k eq 1 select curvelist[1] else curvelist[2];
	      curves := [[Evaluate(f,T): f in ff]: ff in curves];
  	      if IsSquare(Fk) and [x,y] notin curves then
	        pt := [T,x,y];
                printf "Z(17,%o) :  [T,x,y] = %o\n",k,pt;
                assert exists{dd : dd in Table2 | dd[1] eq k and dd[2] eq pt};
              end if;
	    end for;
	  end for;
	end for;
      end for;
    end if;
  end for;
end for;

print "";
IndentPop();

/*

Magma V2.23-8     Fri Jun  4 2021 09:16:36 on chromium [Seed = 115323284]
Type ? for help.  Type <Ctrl>-D to quit.
> load "congr17-check-paper.m";
Loading "congr17-check-paper.m"
Loading "congr17-formulae.m"

Section 1 : Introduction

Section 2 : Verification via modularity

Section 3 : Verification via genus 2 Jacobians

    Computing the Rosenhain invariants
    Time: 4.570
    Computing the genus 2 curve
    Time: 0.600
    Time: 0.940
    Checking the maps C -> E1 and C -> E2

Section 4 : The modular curve X(17)

    Creating the z-curve
    Proposition 4.1 : The curves C1 and C2
    Time: 1.880
    Time: 0.460
    Lemma 4.4 : The map j : X(17) -> P^1
    Time: 0.150
    Time: 1.410
    Time: 0.000

Section 5 : Computations in the symplectic case

    Some random Fp-points on X(17)
    The bi-invariants of degree (2,2)
    The bi-invariants of degree (3,3)
    Creating the surface Sigma_1
    Time: 3.470
    Checking the relations
    Method = modp
    Method = numerical
    Proposition 5.3 : Showing Sigma_1 is birational to Sigma
    Computing the inverse
    Time: 0.170
    Time: 5.830
    Extra relations for computing the j-maps
    Method = modp
    Method = numerical
    Solving for j1 + j2 and j1*j2
    Time: 4.700
    Alternative way to represent the function field Q(Sigma)
    Time: 5.200
    Time: 1.280
    Checking j1 + j2 and j1*j2 against the stored answers
    Time: 0.300
    Time: 0.490

Section 6 : Computations in the anti-symplectic case

    The bi-invariants of degree (2,2)
    The bi-invariants of degrees (3,1) and (1,3)
    Time: 1.480
    The bi-invariants of degree (3,3)
    Creating the surface Sigma_3
    Time: 0.250
    Checking the relations
    Method = modp
    Method = numerical
    Proposition 6.2 : Showing Sigma_3 is birational to Sigma
    Computing the inverse
    Time: 0.040
    Time: 0.060
    Extra relations for computing the j-maps
    The bi-invariants of degrees (3,2) and (2,3)
    Method = modp
    Method = numerical
    Solving for j1 + j2 and j1*j2
    Time: 7.240
    Alternative way to represent the function field Q(Sigma)
    Time: 7.680
    Time: 0.990
    Checking j1 + j2 and j1*j2 against the stored answers
    Time: 0.690
    Time: 0.890

Section 7 : Some modular curves

    Checking Table 1 (including the moduli interpretations)
    m = 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 18 19 20 21 24 25 27 32 36 49 

Section 8 : Examples and further questions

    Some curves on Z(17,1)
    (x,y) = [ 0, -T^3 ]                          at infinity
    (x,y) = [ -T, -T^2 - T ]                     corresponds to X0(18)
    (x,y) = [ -T, -T ]                           spurious
    (x,y) = [ T^2, T^2 ]                         corresponds to X0(25)
    (x,y) = [ T^2, -T^4 + T^2 ]                  spurious
    (x,y) = [ T^3 + T^2, T^4 + 2*T^3 + T^2 ]     at infinity
    (x,y) = [ T^3 - T, T^4 - T^2 - T ]           at infinity
    Some curves on Z(17,3)
    (x,y) = [ -T, -T ]                           at infinity
    (x,y) = [ -T, -T^2 - T ]                     at infinity
    (x,y) = [ T^2, -T^4 + T^2 ]                  at infinity
    (x,y) = [ T^3 + T^2, -T^5 - T^4 + T^2 ]      spurious
    Searching for rational points on Z(17,1) and Z(17,3)
    Z(17,3) :  [T,x,y] = [ 1, 2, 4 ]
    Z(17,1) :  [T,x,y] = [ -3, -27, 108 ]
    Z(17,1) :  [T,x,y] = [ 1/3, -2/75, -11/125 ]
    Z(17,1) :  [T,x,y] = [ -5/6, -5/24, 5/16 ]
    Z(17,3) :  [T,x,y] = [ 9/7, 27/49, -54/49 ]
    Z(17,3) :  [T,x,y] = [ -5/14, 125/392, 375/1568 ]
    Z(17,3) :  [T,x,y] = [ 11/39, 1771/6591, 116380/257049 ]

*/

