//
// This is a MAGMA file accompanying the article :
//
//  "On families of n-congruent elliptic curves" 
//  T. Fisher (May 2011)
// 
// The formulae of the paper are listed (and checked) in the order
// of the paper. Some of the more time consuming checks are commented
// out in this version. The quintic forms in Section 6.5 are loaded in a 
// from a separate file "quintics.m".

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

print "Section 1 : Introduction";

data := [<11,"190b1","2470a1",1>,
         <13,"52a2","988b1",1>,
	 <17,"3675b1","47775b1",-1>];

for d in data do
  n := d[1];
  E1 := EllipticCurve(d[2]);
  E2 := EllipticCurve(d[3]);
  printf "n = %o  E1 = %o  E2 = %o \n",n,d[2],d[3];
  N1 := Conductor(E1);
  N2 := Conductor(E2);
  bp := PrimeDivisors(N1*N2);
  assert forall{p : p in bp | {ReductionType(E,p): E in [E1,E2]} in 
    [{"Good","Split multiplicative"},
     {"Good","Nonsplit multiplicative"},
     {"Nonsplit multiplicative"},
     {"Additive"}]};
  M := LCM(N1,N2);
  mu := M*(&*[1+(1/p) : p in PrimeDivisors(M)]);
  for E in [E1,E2] do
    m,deg := IsogenousCurves(E);
    assert [#m,deg] in {[1,1],[2,2]};
  end for;
  bd := Floor(mu/6) + 1;
  pp := [p : p in [1..bd] | IsPrime(p)];
  v1 := TracesOfFrobenius(E1,bd);
  v2 := TracesOfFrobenius(E2,bd);
  assert forall{i : i in [1..#pp] | pp[i] in bp or (v1[i]-v2[i]) mod n eq 0};
  assert forall{p : p in bp | (N1*N2 mod p^2) eq 0 or
    (v1[Position(pp,p)]*v2[Position(pp,p)]-(p+1)) mod n eq 0};
  assert exists(p){p : p in bp | {ReductionType(E,p): E in [E1,E2]} eq 
    {"Nonsplit multiplicative"}};
  m1 := Valuation(Discriminant(E1),p);
  m2 := Valuation(Discriminant(E2),p); 
  assert LegendreSymbol(m1*m2,n) eq d[4];
end for;

print "Section 1.1 : Some modular curves";

npsl2 := func<n|(n^3/2)*&*[1-1/p^2: p in PrimeDivisors(n)]>;
g := func<n|((n-6)/(12*n))*npsl2(n) + 1>;
print "g(n) =",[g(n): n in [3..17]];

print "Section 1.2 : Statement of results";

K<c4,c6> := FunctionField(Rationals(),2);
R<la,mu> := PolynomialRing(K,2);
DD,cc4,cc6 := HessePolynomials(3,1,[c4,c6]:Variables := [la,mu]);

assert DD eq la^4 - 6*c4*la^2*mu^2 - 8*c6*la*mu^3 - 3*c4^2*mu^4; 
assert cc4 eq c4*la^4 + 4*c6*la^3*mu + 6*c4^2*la^2*mu^2 + 4*c4*c6*la*mu^3 
   + (-3*c4^3 + 4*c6^2)*mu^4;
assert cc6 eq c6*la^6 + 6*c4^2*la^5*mu + 15*c4*c6*la^4*mu^2 
   + 20*c6^2*la^3*mu^3 + 15*c4^2*c6*la^2*mu^4 
   + 6*(3*c4^4 - 2*c4*c6^2)*la*mu^5 + (9*c4^3*c6 - 8*c6^3)*mu^6;

function XE7Equations(a,b:R := PolynomialRing(Universe([a,b]),3))
  x,y,z := Explode([R.i: i in [1..3]]);
  F := a*x^4 + 7*b*x^3*z + 3*x^2*y^2 - 3*a^2*x^2*z^2 - 6*b*x*y*z^2 
    - 5*a*b*x*z^3 + 2*y^3*z + 3*a*y^2*z^2 + 2*a^2*y*z^3 - 4*b^2*z^4;
  G := -a^2*x^4 + 2*a*b*x^3*y - 12*b*x^3*z - (6*a^3 + 36*b^2)*x^2*y^2 
    + 6*a*x^2*z^2 + 2*a^2*b*x*y^3 - 12*a*b*x*y^2*z + 18*b*x*y*z^2 
    + (3*a^4 + 19*a*b^2)*y^4 - (8*a^3 + 42*b^2)*y^3*z + 6*a^2*y^2*z^2 
    - 8*a*y*z^3 + 3*z^4;
  return F,G;
end function;

function XE9Equations(a,b:R := PolynomialRing(Universe([a,b]),4))
  x,y,z,t := Explode([R.i: i in [1..4]]);
  F1 := x^2*t + 6*x*y*z + 6*b*x*t^2 + 6*y^3 - 9*a*y^2*t + 6*a^2*y*t^2 
     - 3*b*z^3 + 3*a^2*z^2*t + 9*a*b*z*t^2  - (a^3 - 12*b^2)*t^3;
  F2 := x^2*z + 6*x*y^2 - 6*a*x*y*t + 2*a^2*x*t^2 - 9*a*y^2*z - 18*b*y*z^2 
     + 12*a^2*y*z*t + a^2*z^3 + 9*a*b*z^2*t - 3*a^3*z*t^2 + a^2*b*t^3;
  G1 := 9*x^2*y + 3*x^2*z - 6*a*x*y*t + 6*b*x*t^2 - 6*a*y^3 + 27*b*y^2*t 
     + 3*a*y*z^2 + 18*b*y*z*t + 3*a^2*y*t^2 + a*z^3 + 3*b*z^2*t 
     + a^2*z*t^2 - a*b*t^3;
  G2 := x^3 + 6*a*x*y*z + 18*b*x*y*t + 3*a*x*z^2 + 6*b*x*z*t + a^2*x*t^2 
     + 9*b*y^3 + 6*a^2*y^2*t - 9*b*y*z^2 + 6*a^2*y*z*t - 3*a*b*y*t^2 
     - 4*b*z^3 + 2*a^2*z^2*t + 2*b^2*t^3;
  return F1,F2,G1,G2;
end function;

function XE11Equations(a,b:R := PolynomialRing(Universe([a,b]),5))
  v,w,x,y,z := Explode([R.i: i in [1..5]]);
  F := v^3 + a*v^2*z - 2*a*v*x^2 + 4*a*v*x*y - 6*b*v*x*z + a*v*y^2 
    + 6*b*v*y*z + a^2*v*z^2 - w^3 + a*w^2*z - 4*a*w*x^2 - 12*b*w*x*z 
    + a^2*w*z^2 - 2*b*x^3 + 3*b*x^2*y + 2*a^2*x^2*z + 6*b*x*y^2 
    + 4*a*b*x*z^2 + b*y^3 - a^2*y^2*z + a*b*y*z^2 + 2*b^2*z^3;
  G := v^2*z + 2*v*w*y + 4*v*x*y + 2*w^2*x - a*w^2*z + 2*w*x^2 
    - 2*a*w*y^2 - 6*b*w*y*z + 6*x^3 - 6*a*x^2*z + 2*a^2*x*z^2 + b*y^3 
    - 2*a^2*y^2*z - 5*a*b*y*z^2 - b^2*z^3;
  return F,G;
end function;

print "Section 1.4 : Introductory example";

K0<z3> := CyclotomicField(3);
K<A,B> := PolynomialRing(K0,2);
R<x,y,z> := PolynomialRing(K,3);
cubic := A*(x^3 + y^3 + z^3) + B*x*y*z;
c4,c6,Delta := Invariants(GenusOneModel(cubic));
D := -27*A^4 - A*B^3;
assert c4 eq -216*A^3*B + B^4;
assert c6 eq 5832*A^6 - 540*A^3*B^3 - B^6;
assert Delta eq D^3;

Z := Integers();
S := Matrix(Z,2,2,[0,1,-1,0]);
T := Matrix(Z,2,2,[1,1,0,1]);
assert (S*T)^3 eq IdentityMatrix(Z,2);

MP := DiagonalMatrix(K,[z3^i: i in [0..2]]);
MQ := Matrix(K,3,3,[<1,3,1>,<2,1,1>,<3,2,1>]);
piS := Matrix(K,3,3,[z3^(i*j): i,j in [0..2]]);
piT := DiagonalMatrix(K,[z3^(-i^2): i in [0..2]]);
assert piS^(-1)*MP*piS eq MQ;
assert piS^(-1)*MQ*piS eq MP^(-1);
assert piT^(-1)*MP*piT eq MP;
assert piT^(-1)*MQ*piT eq z3^2*MP^(-1)*MQ;
assert cubic^piS eq (3*A + B)*(x^3 + y^3 + z^3) + (18*A - 3*B)*x*y*z;
assert cubic^piT eq A*(x^3 + y^3 + z^3) + z3*B*x*y*z;
rhoS := (-1/(z3-z3^2))*Matrix(K0,2,2,[1,1/3,6,-1]);
rhoT := DiagonalMatrix(K0,[z3,z3^2]);
assert rhoS^2 eq -IdentityMatrix(K0,2);
assert (rhoS*rhoT)^3 eq IdentityMatrix(K0,2);
eps := DiagonalMatrix(K0,[1,18]);
assert rhoS^(-1) eq eps*Transpose(rhoS)^(-1)*eps^(-1);
assert rhoT^(-1) eq eps*Transpose(rhoT)^(-1)*eps^(-1);
assert forall{ i : i in [1..3] | 
  f^rhoS eq f and f^rhoT eq z3^i*f where f is [D,c4,c6][i] };

R<la,mu> := PolynomialRing(K,2);
subst1 := [la*A - mu*Derivative(D,B),la*B + mu*Derivative(D,A)];
DD,cc4,cc6 := HessePolynomials(3,1,[c4,c6]:Variables := [la,mu]);
assert cc4 eq Evaluate(c4,subst1);
assert cc6 eq Evaluate(c6,subst1);
subst2 := [la*B + mu*Derivative(D,A),-18*(la*A - mu*Derivative(D,B))];
assert -4*DD*(36*D)^4 eq Evaluate(c4,subst2)*(c4^3-c6^2);
assert -8*cc6*(36*D)^6 eq Evaluate(c6,subst2)*(c4^3-c6^2)^2;

print "Section 2.1 : Klein's equations";

// n = 7
R<a,b,c> := PolynomialRing(Rationals(),3);
aalist := [0,a,b,-c,c,-b,-a];
aa := func<i|aalist[(i mod 7) + 1]>;
M := Matrix(7,7,[aa(i-j)*aa(i+j):i,j in [0..6]]);
F := a^3*b + b^3*c + c^3*a;
assert ideal<R|Pfaffians(M,4)> eq ideal<R|F>;

// n = 9
R<a,b,c,d> := PolynomialRing(Rationals(),4);
aalist := [0,a,-b,d,c,-c,-d,b,-a];
aa := func<i|aalist[(i mod 9) + 1]>;
M := Matrix(9,9,[aa(i-j)*aa(i+j):i,j in [0..8]]);
F1 := a^2*b + b^2*c + c^2*a;
F2 := a*b^2 + b*c^2 + c*a^2 - d^3;
I1 := ideal<R|a,b,c>;
I2 := ideal<R|a - b,b - c,d>;
I3 := ideal<R|a + b + c,b^2 + b*c + c^2,d>;
I4 := ideal<R|F1,F2>;
J := ideal<R|Pfaffians(M,4)>;
assert Radical(J) eq I1 meet I2 meet I3 meet I4;

// n = 11
R<a,b,c,d,e> := PolynomialRing(Rationals(),5);
aalist := [0,a,-c,b,e,d,-d,-e,-b,c,-a];
aa := func<i|aalist[(i mod 11) + 1]>;
M := Matrix(11,11,[aa(i-j)*aa(i+j):i,j in [0..10]]);
F := a^2*b + b^2*c + c^2*d + d^2*e + e^2*a;
H := Matrix(5,5,[Derivative(Derivative(F,i),j): i,j in [1..5]]);
I := ideal<R|Minors(H,4)>;
assert I eq Radical(ideal<R|Pfaffians(M,4)>);
HH := Determinant(H);
assert I eq Radical(ideal<R|[Derivative(HH,i): i in [1..5]]>);

print "Section 3.1 : Formulae in the case n = 7";

K0<z7> := CyclotomicField(7);
g7 := 1 + 2*(z7^3 + z7^5 + z7^6);
M2 := (1/g7)*Matrix(K0,3,3,[z7^(i*j)-z7^(-i*j):i,j in [1,2,4]]);
M3 := Matrix(K0,3,3,[<1,2,1>,<2,3,1>,<3,1,1>]);
M7 := DiagonalMatrix(K0,[z7,z7^4,z7^2]);
G7 := MatrixGroup<3,K0|[M2,M7]>;
assert M3 in G7;
assert #G7 eq 168;

function Level7Invariants(R)
  a,b,c := Explode([R.i: i in [1..3]]);
  F := a^3*b + b^3*c + c^3*a;
  Hmat := Matrix(3,3,[Deriv(Deriv(F,R.i),R.j): i,j in [1..3]]);
  H := (-1/54)*Determinant(Hmat);
  DH := [Deriv(H,R.i):i in [1..3]];
  Tmat := VerticalJoin(Hmat,Matrix(1,3,DH));
  Tmat := HorizontalJoin(Tmat,Matrix(4,1,DH cat [0]));
  c4 := (1/9)*Determinant(Tmat);
  Jmat := Matrix(3,3,[Deriv(f,R.i): i in [1..3],f in [F,H,c4]]);
  c6 := (1/14)*Determinant(Jmat);
  return [F,H,c4,c6];
end function;

R<a,b,c> := PolynomialRing(K0,3);
invariants := Level7Invariants(R);
F,H,c4,c6 := Explode(invariants);
assert F^M2 eq F and F^M7 eq F;
assert c4^3 - c6^2 eq 1728*H^7 - 2048*F^9*H + 256*F^7*c4 + 22016*F^6*H^3 
  - 1088*F^4*H^2*c4 - 60032*F^3*H^5 - 88*F^2*H*c4^2 - 1008*F*H^4*c4;

Nabla := func<f|Vector(3,[Deriv(f,i): i in [1..3]])>;
Cross := func<v,w|Vector(
  [v[2]*w[3]-v[3]*w[2],v[3]*w[1]-v[1]*w[3],v[1]*w[2]-v[2]*w[1]])>;

function Level7Covariants(R)
  F,H,c4,c6 := Explode(Level7Invariants(R));
  dF,dH,dc4,dc6 := Explode([Nabla(f) : f in [F,H,c4,c6]]);
  I22 := 448*F^4*H - 48*F^2*c4 - 2048*F*H^3;
  I24 := 128*F^6 - 160*F^3*H^2 - 236*F*H*c4 - 336*H^4;
  ee1 := I22*Cross(dF,dH) - c4*Cross(dF,dc4) + 12*H^2*Cross(dH,dc4);
  ff1 := I24*Cross(dF,dH) - (16*F^4 - 104*F*H^2)*Cross(dF,dc4) 
                                             + c4*Cross(dH,dc4);
  xx := Vector([R.i: i in [1..3]]);
  ee := Vector([ExactQuotient(f,14*c6): f in Eltseq(ee1)]);
  ff := Vector([ExactQuotient(f,14*c6): f in Eltseq(ff1)]);
  return [xx,ee,ff];
end function;

xx,ee,ff := Explode(Level7Covariants(R));
assert Determinant(Matrix([xx,ee,ff])) eq -c6;
assert Determinant(Matrix([Nabla(F),Nabla(H),Nabla(c4)])) eq 14*c6;

function XE7EquationsRaw(c4,c6:R := PolynomialRing(Parent([c4,c6]),3))
  x,y,z := Explode([R.i: i in [1..3]]);
  Delta := (c4^3 - c6^2)/1728;
  F := 12*x^3*z + 108*x^2*y^2 + 3*c4*x^2*z^2 + 72*c4*x*y^2*z 
    - 108*c4*y^4  - 12*c6*x*y*z^2 + 84*c6*y^3*z + c4^2*x*z^3 
    - 15*c4^2*y^2*z^2 + c4*c6*y*z^3 + 768*Delta*z^4;
  G := 3*x^4 + c4*x^3*z - 18*c4*x^2*y^2 - 3*c6*x^2*y*z + 24*c6*x*y^3 
    + 3*c4^2*x*y^2*z - 9*c4^2*y^4 - c4*c6*y^3*z + 168*Delta*x*z^3 
    + 1728*Delta*y^2*z^2 + 5*c4*Delta*z^4;
  return F,G;
end function;

function TildePolys(F,H,c4,c6:R := PolynomialRing(Parent([F,H,c4,c6]),3))
  x,y,z := Explode([R.i: i in [1..3]]);
  Ftilde := F*x^4 + 12*H^3*x^3*z + (108*H^3 - 6*c4*F)*x^2*y^2 
    - 8*c6*F*x*y^3 + 3*c4*H^3*x^2*z^2 + (72*c4*H^3 + 4128*F^2*H^4 
    + 48*c4*F^3*H - 768*F^5*H^2)*x*y^2*z + (-108*c4*H^3 - 3*c4^2*F 
    - 11376*F^2*H^4 + 32*c4*F^3*H + 3392*F^5*H^2 - 256*F^8)*y^4 
    - 12*c6*H^3*x*y*z^2 + (84*c6*H^3 - 16*c6*F^3*H)*y^3*z 
    + (c4^2*H^3 + 688*F*H^7 + 8*c4*F^2*H^4  - 128*F^4*H^5)*x*z^3 
    + (-15*c4^2*H^3 - 10512*F*H^7  - 384*c4*F^2*H^4   
    + 6144*F^4*H^5 + 96*c4*F^5*H^2 - 768*F^7*H^3)*y^2*z^2 
    + (c4*c6*H^3 - 8*c6*F^2*H^4)*y*z^3 + (768*H^10 - 36*c4*F*H^7 
    - c4^2*F^2*H^4 + 176*F^3*H^8 + 16*c4*F^4*H^5 - 64*F^6*H^6)*z^4;
  Gtilde := (3*H^2 + 28*F^3)*x^4 + (c4*H^2 + 168*F^2*H^3)*x^3*z 
    + (-18*c4*H^2 - 816*F^2*H^3 - 24*c4*F^3 + 192*F^5*H)*x^2*y^2 
    - 3*c6*H^2*x^2*y*z + 24*c6*H^2*x*y^3 + (222*F*H^6 + 24*F^4*H^4)*x^2*z^2 
    + (3*c4^2*H^2 + 3744*F*H^6 - 576*F^4*H^4)*x*y^2*z + (-9*c4^2*H^2 
    - 5184*F*H^6 - 240*c4*F^2*H^3 - 4*c4^2*F^3 + 2240*F^4*H^4 + 64*c4*F^5*H 
    - 256*F^7*H^2)*y^4 + (-c4*c6*H^2 + 8*c6*F^2*H^3)*y^3*z 
    + (168*H^9 + 3*c4*F*H^6 + 24*F^3*H^7)*x*z^3 + (1728*H^9 
    - 78*c4*F*H^6 + 816*F^3*H^7 + 24*c4*F^4*H^4 - 192*F^6*H^5)*y^2*z^2 
    + c6*F*H^6*y*z^3 + (5*c4*H^9 + 35*F^2*H^10 - 4*F^5*H^8)*z^4;
  return Ftilde,Gtilde;
end function;

RR<x,y,z> := PolynomialRing(R,3);
Ftilde,Gtilde := TildePolys(F,H,c4,c6:R:=RR);

covs := [xx,Cross(Nabla(F),Nabla(H)),H*ee];
assert Determinant(Matrix(covs)) eq 72*H^4 - 4*c4*F*H;
subst1 := [&+[covs[i][j]*RR.i: i in [1..3]]: j in [1..3]];
assert Ftilde eq Evaluate(F,subst1);

cons := [Nabla(F),Cross(xx,ee),H^2*Nabla(H)];
assert Determinant(Matrix(cons)) eq 72*H^5 - 4*c4*F*H^2;
subst2 := [&+[cons[i][j]*RR.i: i in [1..3]]: j in [1..3]];
assert Gtilde eq Evaluate(F,subst2);

R<H,c4,c6> := PolynomialRing(Rationals(),3);
RR<x,y,z> := PolynomialRing(R,3);
Ftilde,Gtilde := TildePolys(0,H,c4,c6:R:=RR);
FF,GG := XE7EquationsRaw(c4,c6:R:=RR);
I := ideal<R|c4^3 - c6^2 - 1728*H^7>;
assert Coefficients(Ftilde - H^3*FF) subset I;
assert Coefficients(Gtilde - H^2*GG) subset I;
a := -27*c4;
b := -54*c6;
F,G := XE7Equations(a,b:R:=RR);
assert F eq (1/4)*Evaluate(FF,[6*c4*z-(1/3)*y,x,-18*z]);
assert G eq Evaluate(GG,[9*c4*y+z,3*x,108*y]);
invF := (1/(4*H^3))^3*(72*6*H^4)^4;
invG := (1/H^2)^3*(72*324*H^5)^4;
assert R!invF - (-4*(4*a^3+27*b^2)) in I;
assert R!invG - (16*(4*a^3+27*b^2)^2) in I;

print "Section 3.2 : Formulae in the case n = 9";

K0<z9> := CyclotomicField(9);
M2 := (1/3)*Matrix(K0,4,4,[z9^(i*j)-z9^(-i*j):i,j in [1,7,4,3]]);
M3 := Matrix(K0,4,4,[<1,2,1>,<2,3,1>,<3,1,1>,<4,4,1>]);
M9 := DiagonalMatrix(K0,[z9^i: i in [1,4,7,6]]);
G9 := MatrixGroup<4,K0|[M2,M9]>;
assert M3 in G9;
assert #G9 eq 648;

Nabla := func<f|Vector(4,[Deriv(f,i): i in [1..4]])>;
S4 := SymmetricGroup(4);
MySign := func<i,j,k,l|#{i,j,k,l} lt 4 select 0 else Sign(S4![i,j,k,l])>;
Star := func<mat,vec|Vector([&+[mat[i,j]*vec[k] 
     : i,j,k in [1..4]| MySign(i,j,k,l) eq 1]: l in [1..4]])>;

function Level9Covariants(R)
  a,b,c,d := Explode([R.i: i in [1..4]]);
  _<t> := PolynomialRing(R);
  D := -(a^3 + b^3 + c^3 - 3*a*b*c)*d;
  I6 := 2*(a^5*b + b^5*c + c^5*a) + 5*(a^4*c^2 + b^4*a^2 + c^4*b^2) 
       + 20*a*b*c*(a^2*b + b^2*c + c^2*a) - 3*d^6;
  A := Matrix(4,4,[Deriv(Deriv(D,i),j): i,j in [1..4]]);
  B := Matrix(4,4,[Deriv(Deriv(I6,i),j): i,j in [1..4]]);
  poly := Determinant(A + t*B);
  I12 := (1/2700)*Coefficient(poly,2);
  assert poly eq 81*D^2 - 1620*D*I6*t + 2700*I12*t^2 + 90000*D^4*t^4;
  Ev := Evaluate;
  Lambda := func<q,r|Matrix(4,4,[
    [0,Ev(q,[a,b,c,0])*d^2,-Ev(q,[c,a,b,0])*d^2,Ev(r,[a,b,c,0])],
    [-Ev(q,[a,b,c,0])*d^2,0,Ev(q,[b,c,a,0])*d^2,Ev(r,[b,c,a,0])],
    [Ev(q,[c,a,b,0])*d^2,-Ev(q,[b,c,a,0])*d^2,0,Ev(r,[c,a,b,0])],
    [-Ev(r,[a,b,c,0]),-Ev(r,[b,c,a,0]),-Ev(r,[c,a,b,0]),0]])>;
  Lambda0 := Lambda(3*(a^2 + 2*b*c),2*(a^3*b + c^3*b) - 3*a^2*c^2 - b^4);
  Lambda1 := Lambda(3*(c^2 + 2*a*b),2*(b^3*c + a^3*c) - 3*b^2*a^2 - c^4);
  Lambda2 := Lambda(3*(b^2 + 2*c*a),2*(c^3*a + b^3*a) - 3*c^2*b^2 - a^4);
// We take transposes below to compensate for the fact vectors in Magma 
// are row vectors (not column vectors as in the paper)
  x1 := Vector([a,b,c,d]);
  x7 := Nabla(D)*Transpose(Lambda0);
  e7 := Nabla(D)*Transpose(Lambda1);
  e9 := Nabla(I6)*Transpose(Lambda1);
  f7 := Nabla(D)*Transpose(Lambda2);
  u5 := Star(Lambda1,x1);
  v11 := Star(Lambda2,x7);
  c4 := &+[x7[i]*u5[i]: i in [1..4]];
  c6 := (1/2)*&+[x7[i]*Nabla(I12)[i]: i in [1..4]]; 
  f1 := a^8 + 35*a^6*b*c + 42*a^5*b^3 + 7*a^5*c^3 + 105*a^4*b^2*c^2 
     + 28*a^4*b*d^3 + 231*a^3*b^4*c - 196*a^3*b*c^4 + 14*a^3*c^2*d^3 
     + 14*a^2*b^6 - 70*a^2*b^3*c^3 - 84*a^2*b^2*c*d^3 - 21*a^2*c^6 
     + 28*a*b^4*d^3 - 105*a*b^2*c^5 - 14*a*b*c^3*d^3 - 27*a*c*d^6 
     + 19*b^7*c - 35*b^4*c^4 + 14*b^3*c^2*d^3 + 27*b^2*d^6 - 27*b*c^7 
     + 14*c^5*d^3;
  f2 := Evaluate(f1,[b,c,a,d]);
  f3 := Evaluate(f1,[c,a,b,d]);
  v9 := Vector([f1*d,f2*d,f3*d,-a*f1-b*f2-c*f3]);
  x13 := v9*Transpose(Lambda1);
  x15 := v11*Transpose(Lambda1);
  f13 := v9*Transpose(Lambda0);
  return [D,I6,I12,c4,c6],[x1,x7,x13,x15,e7,e9,f7,f13,u5,v9,v11];
end function;

R<a,b,c,d> := PolynomialRing(K0,4);
F1 := a^2*b + b^2*c + c^2*a;
F2 := a*b^2 + b*c^2 + c*a^2 - d^3;
I := ideal<R|F1,F2>;
invariants,columns := Level9Covariants(R);
D,I6,I12,c4,c6 := Explode(invariants);
assert forall{I : I in [D,I6] | I^M2 eq I};
// time assert forall{I : I in [I12,c4,c6] | I^M2 eq I};
assert forall{I : I in [D,I6,I12,c6] | I^M9 eq I};
assert c4^M9 eq z9^3*c4;
assert c4 notin I;
assert c6 notin I;
assert c4^3 - c6^2 - 1728*D^9 in I;

M2R := ChangeRing(M2,R);
M9R := ChangeRing(M9,R);
types := [0,0,0,0,1,1,2,2,1,2,2];
for i in [1..#columns] do
  x := columns[i];
  r := types[i];
//  time assert Vector([x[j]^M2: j in [1..4]]) eq x*M2R;
  assert Vector([x[j]^M9: j in [1..4]]) eq z9^(3*r)*x*M9R;
  if i eq 8 then M2R := M2R^(-1); M9R := M9R^(-1); end if;
end for;

x1,x7,x13,x15,e7,e9,f7,f13,u5,v9,v11 := Explode(columns);

// Proof of Lemma 3.12 :

// The covariant columns of degree = 1 mod 4 form a free K[D,c4]-module
// of rank 4 with basis
aa := [x1,e9,f13,x13];
// The covariant columns of degree = 3 mod 4 form a free K[D,c4]-module
// of rank 4 with basis
bb := [x7,e7,f7,x15];
// The contravariant columns of degree = 1 mod 4 form a free K[D,c4]-module
// of rank 4 with basis
cc := [Nabla(I6),u5,v9,Nabla(c6)];
// The contravariant columns of degree = 3 mod 4 form a free K[D,c4]-module
// of rank 4 with basis
dd := [Nabla(D),Nabla(c4),v11,Nabla(I12)];

// The above claims are proved along the same lines as Lemma 3.5 
// by considering the following dot products
M := Matrix(4,4,[0,0,0,-2*c6,0,0,-c6,0,0,2*c6,0,0,18*c6,0,0,108*D^3*c6]);
N := Matrix(4,4,[0,0,0,c6,0,0,-12*c6,0,0,-c6,0,0,2*c6,0,0,-6*D^2*c6]);
for i,j in [1..4] do
  assert (&+Eltseq(cc[i]*aa[j]) - M[i,j]) in I;
  assert (&+Eltseq(dd[i]*bb[j]) - N[i,j]) in I;
end for;

// Hilbert series calculations showing that the bases in Lemma 3.12 
// consist of elements of the right degree.
P<t> := FunctionField(Rationals());
// Covariant columns :
f1 := (t + t^9 + t^13 + t^13)/((1-t^6)*(1-t^12));   // deg = 1 mod 4
f3 := (t^7 + t^7 + t^7 + t^15)/((1-t^6)*(1-t^12));  // deg = 3 mod 4
g0 := (t + t^7 + t^13 + t^15)/((1-t^6)*(1-t^18));    // type 0
g1 := (t^7 + t^9 + t^13 + t^19)/((1-t^6)*(1-t^18));  // type 1
g2 := (t^7 + t^13 + t^19 + t^21)/((1-t^6)*(1-t^18)); // type 2
assert f1 + f3 eq g0 + g1 + g2;
// Contravariant columns :
f1 := (t^5 + t^5 + t^9 + t^17)/((1-t^6)*(1-t^12));  // deg = 1 mod 4
f3 := (t^3 + t^11 + t^11 + t^11)/((1-t^6)*(1-t^12));// deg = 3 mod 4
g0 := (t^3 + t^5 + t^11 + t^17)/((1-t^6)*(1-t^18));  // type 0
g1 := (t^5 + t^11 + t^15 + t^17)/((1-t^6)*(1-t^18)); // type 1
g2 := (t^9 + t^11 + t^17 + t^23)/((1-t^6)*(1-t^18)); // type 2
assert f1 + f3 eq g0 + g1 + g2;

// Checking the bases in Lemma 3.12 are independent
assert Determinant(Matrix([x1,x7,x13,x15])) + c4^3 in I;
assert Determinant(Matrix([e7,e9,x1,x7])) + 2*c4^2 in I;
assert Determinant(Matrix([f7,f13,e7,e9])) + 2*c4^3 in I;
assert Determinant(Matrix([Nabla(f): f in [D,I6,I12,c6]])) - 72*c4^3 in I;
assert Determinant(Matrix([u5,Nabla(c4),Nabla(D),Nabla(I6)])) - 24*c4^2 in I;
assert Determinant(Matrix([v9,v11,u5,Nabla(c4)])) - 12*c4^3 in I;

function XE9EquationsRaw(c4,c6:R := PolynomialRing(Parent([c4,c6]),4))
  x,y,z,t := Explode([R.i: i in [1..4]]);
  F1 := 24*x^2*t - 96*x*y*z + 64*y^3 + 24*c4*y^2*t + 48*c4*y*z^2 
     - 48*c6*y*z*t + 12*c4^2*y*t^2 - 8*c6*z^3 + 12*c4^2*z^2*t 
     - 6*c4*c6*z*t^2 - (c4^3 - 2*c6^2)*t^3;
  F2 := 16*x^2*z - 64*x*y^2 - 16*c4*x*y*t - 16*c4*x*z^2 + 16*c6*x*z*t 
     - 4*c4^2*x*t^2 + 80*c4*y^2*z - 32*c6*y^2*t - 32*c6*y*z^2 
     + 32*c4^2*y*z*t - 8*c4*c6*y*t^2 + 8*c4^2*z^3 - 12*c4*c6*z^2*t 
     + (2*c4^3 + 4*c6^2)*z*t^2 - c4^2*c6*t^3;
  G1 := -72*x^2*y - 144*x^2*z + 24*c4*x*y*t - 48*c4*x*z*t - 8*c6*x*t^2 
     - c4*y^3 + 18*c4*y^2*z + 3*c6*y^2*t + 180*c4*y*z^2 + 12*c6*y*z*t 
     - 3*c4^2*y*t^2 + 792*c4*z^3 + 12*c6*z^2*t + 2*c4^2*z*t^2 + c4*c6*t^3;
  G2 := -864*x^3 + 216*c4*x^2*t + 2592*c4*x*y*z + 216*c6*x*y*t 
     + 15552*c4*x*z^2 + 432*c6*x*z*t - 72*c4^2*x*t^2 - 9*c6*y^3 
     + 162*c6*y^2*z + 27*c4^2*y^2*t + 4212*c6*y*z^2 + 108*c4^2*y*z*t 
     - 27*c4*c6*y*t^2 + 26136*c6*z^3 + 972*c4^2*z^2*t + 18*c4*c6*z*t^2 
     + (5*c4^3 + 4*c6^2)*t^3;
  return F1,F2,G1,G2;
end function;

A := d^3;
B := a^3 + b^3 + c^3 + 6*a*b*c;
alpha := (1/9)*Matrix(2,2,[3*B,9*A,-(108*A^3+B^3),9*A*B^2]);
beta := Matrix(2,2,[6*A,-B,162*A*B^2,9*(108*A^3 + B^3)]);
assert Determinant(alpha) + (4/9)*D^3 in I;
assert Determinant(beta) + 216*D^3 in I;
RR<x,y,z,t> := PolynomialRing(R,4);
covs := [x1,f7,D*e9,D*x15];
cons := [u5,D^2*Nabla(D),Nabla(I12),D^2*v9];
assert Determinant(Matrix(covs)) - 3456*D^10 in I;
assert Determinant(Matrix(cons)) - 1152*D^11 in I;
subst1 := [&+[covs[i][j]*RR.i: i in [1..4]]: j in [1..4]];
subst2 := [&+[cons[i][j]*RR.i: i in [1..4]]: j in [1..4]];
FF := [Evaluate(F,subst1): F in [F1,F2]];
FF := [&+[alpha[i,j]*FF[j]: j in [1..2]]: i in [1..2]];
GG := [Evaluate(F,subst2): F in [F1,F2]];
GG := [&+[beta[i,j]*GG[j]: j in [1..2]]: i in [1..2]];
F1,F2,G1,G2 := XE9EquationsRaw(c4,c6:R:=RR);
assert Coefficients(FF[1] - D^6*F1) subset I;
assert Coefficients(FF[2] - D^6*F2) subset I;
assert Coefficients(GG[1] - D^6*G1) subset I;
assert Coefficients(GG[2] - D^6*G2) subset I;

R<c4,c6> := PolynomialRing(Rationals(),2);
RR<x,y,z,t> := PolynomialRing(R,4);
FF1,FF2,GG1,GG2 := XE9EquationsRaw(c4,c6:R:=RR);
a := -27*c4;
b := -54*c6;
F1,F2,G1,G2 := XE9Equations(a,b:R:=RR);
subst := [x - 27*c4*z - 162*c6*t,9*y + 81*c4*t,-54*z,324*t];
assert F1 eq (1/7776)*Evaluate(FF1,subst);
assert F2 eq -(1/864)*Evaluate(FF2,subst);
subst := [-x + 9*c4*t,36*y + 9*z,(3/2)*z,108*t];
assert G1 eq -(1/288)*Evaluate(GG1,subst);
assert G2 eq (1/864)*Evaluate(GG2,subst);
Delta := (c4^3-c6^2)/1728;
invF := ((4/9)/(7776*864))^6*(3456*9*(-54)*324)^9*Delta^4;
invG := (216/(288*864))^6*(1152*(-36)*(3/2)*108)^9*Delta^5;
assert invF eq -2^10*(4*a^3 + 27*b^2)^4;
assert invG eq 2^8*(4*a^3 + 27*b^2)^5;

print "Section 3.3 : Formulae in the case n = 11";

K0<z11> := CyclotomicField(11);
g11 := 1 + 2*&+[z11^i: i in [1,3,9,5,4]];
M2 := (1/g11)*Matrix(K0,5,5,[z11^(i*j)-z11^(-i*j):i,j in [1,3,9,5,4]]);
M5 := Matrix(K0,5,5,[<1,2,1>,<2,3,1>,<3,4,1>,<4,5,1>,<5,1,1>]);
M11 := DiagonalMatrix(K0,[z11^i : i in [1,9,4,3,5]]);
G11 := MatrixGroup<5,K0|[M2,M11]>;
assert M5 in G11;
assert #G11 eq 660;

function Level11Invariants(R)
  a,b,c,d,e := Explode([R.i: i in [1..5]]);
  K := BaseRing(R);
  M5 := Matrix(K,5,5,[<1,2,1>,<2,3,1>,<3,4,1>,<4,5,1>,<5,1,1>]);
  CyclicSum := func<f|&+[f^(M5^i): i in [0..4]]>;
  F := CyclicSum(a^2*b);
  H := 3*a*b*c*d*e + CyclicSum(a^3*c^2 - a^3*d*e);
  I7 := CyclicSum(a^6*e + 3*a^5*d^2 - 15*a^4*b*c*e + 5*a^3*b^3*d 
     + 15*a^3*b*c*d^2);
  I8 := CyclicSum(a^7*c - 7*a^4*b*d^3 - 7*a^4*d*e^3 + 7*a^3*b^2*c^3 
     + 21*a^3*c^2*d^2*e);
  I9 := CyclicSum(5*a^8*d + a^6*b^3 - 63*a^6*b*c*d + 40*a^5*c^3*e 
     + 17*a^5*c*d*e^2 - 10*a^4*b^4*c + 163*a^4*b^2*c^2*d 
     + 86*a^4*b^2*d^2*e - 123*a^4*c^3*d^2 - 23*a^4*c*d^3*e 
     - 97*a^3*b^3*c*e^2 - 51*a^3*b*c^2*d*e^2);
  c4 := -CyclicSum(a^11 - 16*a^9*b*c - 6*a^8*c*e^2 + 93*a^7*b^2*c^2 
     - 54*a^7*b^2*d*e - 285*a^7*c^3*d + 243*a^7*c*d^2*e + 261*a^6*b^2*d^3 
     - 6*a^6*b*c^2*e^2 + 153*a^6*b*d*e^3 - 234*a^6*c*d^4 - 30*a^5*b^5*e 
     - 267*a^5*b^3*c^3 + 372*a^5*b^3*c*d*e + 738*a^5*b*c^4*d 
     - 1143*a^5*b*c^2*d^2*e - 678*a^5*b*d^3*e^2 + 9*a^5*c^2*e^4 
     + 225*a^4*b^4*e^3 - 858*a^4*b^3*c*d^3 + 204*a^4*b^2*c^3*e^2 
     - 681*a^4*b^2*c*d*e^3 + 714*a^4*b*c^2*d^4 - 108*a^4*c^2*d^2*e^3 
     + 1224*a^3*b^2*c^3*d^2*e);
  return [F,H,I7,I8,I9,c4];
end function; 

R<a,b,c,d,e> := PolynomialRing(K0,5);
invariants := Level11Invariants(R);
// assert forall{I : I in invariants | I^M2 eq I and I^M11 eq I};
F,H,I7,I8,I9,c4 := Explode(invariants);
A := Matrix(5,5,[Deriv(Deriv(F,i),j): i,j in [1..5]]);
B := Matrix(5,5,[Deriv(Deriv(H,i),j): i,j in [1..5]]);
_<t> := PolynomialRing(R);
poly := Determinant(A+t*B);
truncpoly := &+[Coefficient(poly,i)*t^i: i in [0..3]];
assert truncpoly eq 32*H - 32*I7*t - 24*I9*t^2 - 8*c4*t^3;

c6tilde := a^9*b^10 - 509*b^18*d - 14107*b^14*d^4*e + 510*b^9*c^10 
  + 42326*b^7*d^12 + 20669*b^3*d^15*e - 14107*b^2*d^2*e^15 
  - 277419*b*c^2*d^10*e^6 - 248909*b*c*d^16*e - 209926*b*c*d^5*e^12 
  + 762409*b*d^11*e^7 + b*e^18 - 1018*c^18*e - 14107*c^16*d*e^2 
  - 586835*c^12*d^3*e^4 + 197780*c^10*d^4*e^5 + 1019*c^9*d^10 
  - 787130*c^8*d^5*e^6 + 15634*c^7*d^11*e + 42326*c^7*e^12 
  + 2007576*c^6*d^6*e^7 + 247382*c^5*d^12*e^2 - 528424*c^5*d*e^13 
  - 616653*c^4*d^7*e^8 + 376744*c^3*d^13*e^3 + 1067732*c^3*d^2*e^14 
  - 225004*c^2*d^8*e^9 + 463659*c*d^14*e^4 - 582142*c*d^3*e^15 
  + 70511*d^9*e^10;

Hmat := Matrix(5,5,[Derivative(Derivative(F,i),j): i,j in [1..5]]);
I := ideal<R|Minors(Hmat,4)>;
assert c6tilde^2 - a*b*c*d*e*(c4^3 - 1728*F^11) in I;

Nabla := func<f|Vector(5,[Deriv(f,i): i in [1..5]])>;

function Alt2Mat(M)
  n := Nrows(M);
  ii := [[i,j]: i,j in [1..n] | i lt j];
  MM := Matrix(#ii,#ii,[[MM[j[1],j[2]] : j in ii] 
      where MM is Transpose(M)*E*M 
      where E is Matrix(n,n,[<k[1],k[2],1>,<k[2],k[1],-1>]): k in ii]);
  return MM;
end function; 

function Level11Covariants(R)
  F,H,I7,I8,I9,c4 := Explode(Level11Invariants(R));
  a,b,c,d,e := Explode([R.i: i in [1..5]]);
  f4 := 2*a^2*e^2 + 4*a*b^2*c - 4*a*c^2*d + 4*b*c*e^2 + d^4;
  f5 := -5*a^3*c*e + 5*a^2*b^2*d + 5*a^2*c*d^2 + 5*a*b*c^2*e 
     - 10*a*b*d*e^2 + b^5 - 5*b^3*c*d + 5*b*d^3*e + 5*c^2*e^3;
  f9 := -14*a^6*b*d*e - 8*a^5*b*d^3 + 9*a^5*c^2*e^2 + 2*a^5*d*e^3 
     + 8*a^4*b^4*e + 5*a^4*b^2*c^3 + 63*a^4*b^2*c*d*e + 6*a^4*c^4*d 
     - 18*a^4*c^2*d^2*e + 8*a^4*d^3*e^2 + 31*a^3*b^4*d^2 
     - 21*a^3*b^3*e^3 + 47*a^3*b^2*c*d^3 + 35*a^3*b*c^3*e^2 
     + 14*a^3*b*c*d*e^3 - 12*a^3*c^2*d^4 + 10*a^3*d^5*e + 3*a^2*b^5*c*e 
     - 26*a^2*b^3*c^4 - 42*a^2*b^3*c^2*d*e - 75*a^2*b^3*d^2*e^2 
     + 3*a^2*b^2*e^5 + 18*a^2*b*c^5*d - 30*a^2*b*c^3*d^2*e 
     - 36*a^2*b*c*d^3*e^2 + 2*a^2*c^3*e^4 - 9*a^2*c*d*e^5 + a^2*d^7 
     - 2*a*b^7*d - 6*a*b^5*c*d^2 + 50*a*b^4*c*e^3 - 7*a*b^3*c^2*d^3 
     - 6*a*b^3*d^4*e - 54*a*b^2*c^4*e^2 - 3*a*b^2*c^2*d*e^3 
     - 9*a*b^2*d^2*e^4 - 29*a*b*c^3*d^4 + 21*a*b*c*d^5*e + a*b*e^7 
     + 9*a*c^5*d*e^2 + 25*a*c^3*d^2*e^3 - 7*a*c*d^3*e^4 - 10*b^6*c^2*e 
     - 2*b^6*d*e^2 + 4*b^4*c^5 + 40*b^4*c^3*d*e - 6*b^4*c*d^2*e^2 
     + 13*b^3*c*e^5 - 3*b^3*d^6 - 15*b^2*c^4*d^2*e - 54*b^2*c^2*d^3*e^2 
     + 31*b^2*d^4*e^3 - 11*b*c^4*e^4 + 3*b*c^2*d*e^5 - 2*b*c*d^7 
     - 7*b*d^2*e^6 - c^7*d^2 + 5*c^5*d^3*e - 9*c^3*d^4*e^2 
     + 8*c*d^5*e^3 - e^9;
  K := BaseRing(R);
  M5 := Matrix(K,5,5,[<1,2,1>,<2,3,1>,<3,4,1>,<4,5,1>,<5,1,1>]);
  x1 := Vector([a,b,c,d,e]);
  x4 := Vector([f4^(M5^i): i in [0..4]]);
  x5 := Vector([f5^(M5^i): i in [0..4]]);
  x9 := Vector([f9^(M5^i): i in [0..4]]);
  dF := Nabla(F);
  dI7 := Nabla(I7);
  Xi0 := Matrix(5,5,[dF[i]*dI7[j] - dF[j]*dI7[i]: i,j in [1..5]]);
  r := func<i,j|(((i-j-2)^3 + i + 2) mod 5 + 1)>;
  s := func<i,j|(((j-i-2)^3 + j + 2) mod 5 + 1)>;
  Xi := Matrix(5,5,[Xi0[r(i,j),s(i,j)] : i,j in [1..5]]);
  x14 := Nabla(I7)*Transpose(Xi);
// Explaining where r and s come from:
  IdMat := IdentityMatrix(Rationals(),10);
  cob := IdMat;
  ii := [[i,j]: i,j in [1..5]|i lt j];
  for k in [1..10] do
    u := ii[k];
    v := [r(u[1],u[2]),s(u[1],u[2])];
    l := Position(ii,v);
    if l ne 0 then 
      cob[k] := IdMat[l];
    else
      l := Position(ii,[v[2],v[1]]);
      cob[k] := -IdMat[l];
    end if;
  end for;
  assert forall{M : M in [M2,M5,M11] | 
    cob*Alt2Mat(M) eq Alt2Mat(Transpose(M^(-1)))*cob};
  return [x1,x4,x5,x9,x14];
end function;

covs := Level11Covariants(R);
M2R := ChangeRing(M2,R);
M11R := ChangeRing(M11,R);
for i in [1..3] do 
  x := covs[i];
  assert Vector([x[j]^M2: j in [1..5]]) eq x*M2R;
  assert Vector([x[j]^M11: j in [1..5]]) eq x*M11R;
end for;

function XE11EquationsRaw(c4,Delta:R := PolynomialRing(Parent([c4,Delta]),5))
  v,w,x,y,z := Explode([R.i: i in [1..5]]);
  F := v^3 + 3*v^2*w + c4*v^2*y + 3*v*w^2 + 2*c4*v*w*y - c4*Delta*v*x^2 
    + 48*Delta*v*x*y + 9*w^3 + 5*c4*w^2*y - c4^2*w^2*z + c4^2*w*y^2 
    - 576*Delta*w*y*z + 72*c4*Delta*w*z^2 - 4*Delta^2*x^3 
    - 72*Delta^2*x^2*z + 4*c4*Delta*x*y^2 + 2*c4^2*Delta*x*y*z 
    - (c4^3*Delta - 1728*Delta^2)*x*z^2 + 64*Delta*y^3 
    - 72*c4*Delta*y^2*z + 12*c4^2*Delta*y*z^2 
    + (c4^3*Delta - 3456*Delta^2)*z^3;
  G := 5*v^3 - c4*v^2*x - 60*v^2*y + 28*c4*v^2*z - 2*c4*Delta*v*w^2 
    - 48*Delta*v*w*x - 240*Delta*v*w*z - 16*c4*v*x*y + 1680*v*y^2 
    - 872*c4*v*y*z + 121*c4^2*v*z^2 + 8*Delta^2*w^3 + 44*c4*Delta*w^2*y 
    - 11*c4^2*Delta*w^2*z + c4*Delta*w*x^2 + 336*Delta*w*x*y 
    - 122*c4*Delta*w*x*z + 25*c4^2*w*y^2 - 14160*Delta*w*y*z 
    + 817*c4*Delta*w*z^2 - 20*Delta*x^3 + 5*c4^2*x^2*y - 1884*Delta*x^2*z
    - 364*c4*x*y^2 + 160*c4^2*x*y*z - 34764*Delta*x*z^2 + 19840*y^3 
    - 10268*c4*y^2*z + 1643*c4^2*y*z^2 - 129220*Delta*z^3;
  return F,G;
end function;

cons := [Nabla(f): f in [F,I7,I8,I9,c4]];
assert Determinant(Matrix(covs)) - (c4^3 - 1728*F^11) in I;
assert Determinant(Matrix(cons)) - 55*(c4^3 - 1728*F^11) in I;

function TildePolys(F,c4:R := PolynomialRing(Parent([F,c4]),5))
  v,w,x,y,z := Explode([R.i: i in [1..5]]);
  Ftilde := F*v^3 + 3*F^2*v^2*w + c4*v^2*y + 3*F^3*v*w^2 + 2*F*c4*v*w*y 
    - c4*v*x^2 + 48*F^5*v*x*y + 9*F^4*w^3 + 5*F^2*c4*w^2*y - c4^2*w^2*z 
    + c4^2*w*y^2 - 576*F^9*w*y*z + 72*F^7*c4*w*z^2 - 4*F^5*x^3 
    - 72*F^8*x^2*z + 4*F^4*c4*x*y^2 + 2*F^2*c4^2*x*y*z 
    + (1728*F^11 - c4^3)*x*z^2 + 64*F^9*y^3 - 72*F^7*c4*y^2*z 
    + 12*F^5*c4^2*y*z^2 + (-3456*F^14 + F^3*c4^3)*z^3;
  Gtilde := 5*F^2*v^3 - c4*v^2*x - 60*F^4*v^2*y + 28*F*c4*v^2*z 
    - 2*F*c4*v*w^2 - 48*F^5*v*w*x - 240*F^6*v*w*z - 16*F^2*c4*v*x*y 
    + 1680*F^6*v*y^2 - 872*F^3*c4*v*y*z + 121*c4^2*v*z^2 + 8*F^6*w^3 
    + 44*F^3*c4*w^2*y - 11*c4^2*w^2*z + F^3*c4*w*x^2 + 336*F^7*w*x*y 
    - 122*F^4*c4*w*x*z + 25*c4^2*w*y^2 - 14160*F^8*w*y*z + 817*F^5*c4*w*z^2 
    - 20*F^7*x^3 + 5*c4^2*x^2*y - 1884*F^8*x^2*z - 364*F^4*c4*x*y^2 
    + 160*F*c4^2*x*y*z - 34764*F^9*x*z^2 + 19840*F^8*y^3 
    - 10268*F^5*c4*y^2*z + 1643*F^2*c4^2*y*z^2 - 129220*F^10*z^3;
  return Ftilde,Gtilde;
end function;

RR<v,w,x,y,z> := PolynomialRing(R,5);
Ftilde,Gtilde := TildePolys(F,c4:R:=RR);

subst1 := [&+[covs[i][j]*RR.i: i in [1..5]]: j in [1..5]];
subst2 := [&+[cons[i][j]*RR.i: i in [1..5]]: j in [1..5]];

print "Checking cubic form for X_E(11)";
assert Coefficients(Evaluate(F,subst1) - Ftilde) subset I;

print "Checking cubic form for X_E^-(11)";
assert Coefficients(Evaluate(F,subst2) - Gtilde) subset I;

R<F,c4,c6> := PolynomialRing(Rationals(),3);
RR<v,w,x,y,z> := PolynomialRing(R,5);
Ftilde,Gtilde := TildePolys(F,c4:R:=RR);
FF,GG := XE11EquationsRaw(c4,F^11:R:=RR);
assert F^4*FF eq Evaluate(Ftilde,[F*v,w,F^7*x,F^2*y,F^4*z]);
assert F^8*GG eq Evaluate(Gtilde,[F^2*v,F^8*w,F^4*x,y,F^3*z]);

FF,GG := XE11EquationsRaw(c4,(c4^3-c6^2)/1728:R:=RR);
v1 := c6*v + 2*c6*w - 6*c4^2*x + 3*c4^2*y - 9*c4*c6*z;
w1 := c6*v + 6*c4^2*x + 3*c4^2*y + 9*c4*c6*z;
v11 := 44*(2*c4*v - 6*c6*w + 33*c6*x + 135*c4^2*y + 810*c4*c6*z);
x11 := 60*(5*v + 729*c4*y + 2187*c6*z);
y11 := 11*(c4*v - 3*c6*w - 6*c6*x);
z11 := 60*(v + 27*c4*y + 81*c6*z);
a := -27*c4;
b := -54*c6;
F,G := XE11Equations(a,b:R:=RR);
subst1 := [-v1,w1,-864*x,-36*c4*x - 108*c6*z,72*y];
subst2 := [v11,-427680*y,x11,-y11,-z11];
det1 := Determinant(Matrix(5,5,[MC(s,RR.i): i in [1..5],s in subst1]));
det2 := Determinant(Matrix(5,5,[MC(s,RR.i): i in [1..5],s in subst2]));
assert Evaluate(FF,subst1) eq 8*c6^3*F;
assert Evaluate(GG,subst2) eq 2^5*3^6*5^3*11^3*c6^3*G;
Delta := (c4^3-c6^2)/1728;
invF := (det1*c6^2)^3/(8*c6^3)^5*Delta^2;
invG := (det2*55*c6^2)^3/(2^5*3^6*5^3*11^3*c6^3)^5*Delta;
assert invF eq -4*(4*a^3 + 27*b^2)^2;
assert invG eq 8*(4*a^3 + 27*b^2);

print "Section 4 : Diagonal Twists";

X := EllipticCurve([0,-1,1,0,0]);
K<la,nu> := FunctionField(X);
assert nu^2 + nu eq la^3 - la^2;

D7 := func<la|
  EllipticCurve([-(la^2-la-1),-(la^3-la^2),-(la^3-la^2),0,0])>;
D9 := func<la|
  EllipticCurve([la^3-3*la^2+4*la-1,la*(la-1)*(la^2-la+1),
                       la*(la-1)^4*(la^2-la+1),0,0])>;
D11 := func<la,nu|
  EllipticCurve([la*nu+2*la-(nu+1)^2,-la*nu*(nu+1)*(la-nu-1),
		      -la^2*nu*(nu+1)*(la-nu-1),0,0])>;

function DD(n:lala := n in [7,9] select [la] else [la,nu])
  case n :
    when 7  : return D7(lala[1]);
    when 9  : return D9(lala[1]);
    when 11 : return D11(lala[1],lala[2]);
  end case;
end function;

assert forall{n : n in [7,9,11] | Order(DD(n)![0,0]) eq n};

function CC(n:lala := n in [7,9] select [la] else [la,nu])
  D := DD(n:lala := lala);
  P := D![0,0];
  a1,a2,a3,a4,a6 := Explode(aInvariants(D));
  b2,b4,b6 := Explode(bInvariants(D));
  m := Integers()!((n-1)/2);
  s0 := m;
  s1 := &+[(i*P)[1]: i in [1..m]];
  s2 := &+[(i*P)[1]^2: i in [1..m]];
  s3 := &+[(i*P)[1]^3: i in [1..m]];
  t := 6*s2 + (a1^2 + 4*a2)*s1 + a1*a3*s0;
  w := 10*s3 + 2*(a1^2 + 4*a2)*s2 + 3*a1*a3*s1 + a3^2*s0;
  AA := [a1,a2,a3,-5*t,-b2*t-7*w];
  return EllipticCurve(AA);
end function;

MyDisc := func<n|case<n| 
  7  : la*(la-1)*(la^3-8*la^2+5*la+1)^7, 
  9  : la*(la-1)*(la^2-la+1)^3*(la^3-6*la^2+3*la+1)^9,
  11 : la*(la-1)*(nu+1)^6*(la*nu + 2*la^2 - 2*la + 1)*f^11 
    where f is ((-3*la + 2)*nu - la^3 + 5*la^2 - 5*la + 1)/(la-1), 
  default : 0 >>;

assert forall{ n : n in [7,9,11] | Discriminant(CC(n)) eq MyDisc(n) };

print "Computing q(lambda)";

q := func<n|case<n| 
  7  : la^4*(la-1), 
  9  : la*(la-1)^7*(la^2-la+1)^3,
  11 : la*nu^2*(la-1)*(la-nu-1)^3, 
  default : 0 >>;

for n in [7,9,11] do
  D := DD(n);
  qq := q(n);
  P := D![0,0];
  flag,f := IsPrincipal(n*(Divisor(P)-Divisor(D!0)));
  assert flag;
  offset := Evaluate(f,-P)/qq;
  r := Integers()!((n-1)/2);
  for m in [2..n-1] do
    r := Integers()!((n-1)/2);
    m1 := ((m + r) mod n) - r;
    assert IsPower(Evaluate(f,m*P)*(qq^m1)/offset,n);
  end for;
end for;

// Remark 4.5 ( = alternative proof of Corollary 4.4 for n = 9)

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

A0 := 1/(2^12*3^14*(27*la-1)^5)*Matrix(K,2,2,[ -27, -3, 324*la - 3, 1 ]);
B0 := Matrix(K,4,4,[ 472392*la^2 + 43740*la - 81, 52488*la^2 - 3888*la - 9,
  -324*la + 3, 1, -2834352*la^2 - 262440*la + 486, 8748*la + 162, 
  972*la + 18, -6, 17006112*la^2 - 551124*la + 16767, -26244*la + 243, 
  -81, 9, -708588*la + 6561, 26244*la - 243, 81, -9 ]);
kk := (la^2-la+1)^5*(la^3 + 3*la^2 - 6*la + 1)*(la^3 - 6*la^2 + 3*la + 1)^3;
A1 := (1/kk)*Matrix(K,2,2,[ la^2 - la + 1, 0, 0, -la^3 - 3*la^2 + 6*la - 1 ]);
M := Matrix(K,3,3,[ la - 2, 2*la - 1, la^2 - 1, la + 1, la^2 - 2*la, 
      2*la^2 - la, 2*la - 1, la^2 - 1, la^3 - 3*la^2 + 2*la ]);
M := DiagonalJoin(M,IdentityMatrix(K,1));
B1 := M*DiagonalMatrix(K,[(la^3+3*la^2-6*la+1)^2,
   (la^3+3*la^2-6*la+1)*(la^2-la+1),(la^2-la+1)^2,
   -(la^3-6*la^2+3*la+1)*(la^2-la+1)^2]);
kk := 2^16*3^12*(27*la-1)^5;
A2 := (1/kk)*Matrix(K,2,2,[ 972*la + 45, 9, -810*la + 3, -54*la - 1 ]);
B2 := Matrix(K,4,4,[ 2916, -324, 648, -36, -209952*la^2 + 3888*la - 18, 
  648*la - 6, -1296*la + 12, 2, -17496*la + 162, 1944*la - 18, 
  -1944*la - 36, 6, -52488*la + 486, 1944*la + 90, -5832*la - 108, 18 ]);
kk := (la^2 - la + 1)^5*(la^3 - 6*la^2 + 3*la + 1)^3*(la^3 + 3*la^2 - 6*la + 1)^2;
A3 := (1/kk)*Matrix(K,2,2,[ -9*la^8 + 36*la^7 - 81*la^6 + 117*la^5 - 117*la^4 
  + 81*la^3 - 36*la^2 + 9*la, -la^9 - 9*la^8 - 9*la^7 + 78*la^6 + 36*la^5 
  - 315*la^4 + 321*la^3 - 117*la^2 + 18*la - 1, -la^8 - 5*la^7 + 8*la^6 
  + 25*la^5 - 73*la^4 + 88*la^3 - 55*la^2 + 13*la - 1, -3*la^8 - 15*la^7
  + 24*la^6 + 75*la^5 - 219*la^4 + 264*la^3 - 165*la^2 + 39*la - 3 ]);
M := Matrix(K,3,3,[ -la^2 + la, 1, la, la - 1, 1, -la + 1, la, -1, 1 ]);
M := DiagonalJoin(M,IdentityMatrix(K,1));
B3 := M*DiagonalMatrix(K,[(la^2-la+1)^2,(la^3+3*la^2-6*la+1)^2,
   (la^2-la+1)*(la^3+3*la^2-6*la+1), -(la^2-la+1)*(la^3+3*la^2-6*la+1) ]);

function Apply(tr,polys)
  polys := [f^Transpose(tr[2]): f in polys];
  polys := [&+[tr[1][i,j]*polys[j]: j in [1,2]]: i in [1,2]];
  return polys;
end function;

C3 := EllipticCurve([1,0,la,-5*la,-la*(7*la + 1)]); // mu3 - family
C9 := CC(9:lala := [la]);                           // mu9 - family 
la3 := la*(la-1)*(la^2-la+1)^3/(la^3 + 3*la^2 - 6*la + 1)^3;
c4,c6 := Explode(cInvariants(C3));
a := -27*c4;
b := -54*c6;
u := la^3 + 3*la^2 - 6*la + 1;
assert cInvariants(C9) eq [u^4*Evaluate(c4,la3),u^6*Evaluate(c6,la3)];

R<x,y,z,t> := PolynomialRing(K,4);
F1,F2,G1,G2 := XE9Equations(a,b:R:=R);
mons := MD(R,3);

polys := Apply(<A0,B0>,[F1,F2]);
assert polys eq [la*x^2*y + y^2*z + z^2*x - 3*(z^2*y + la*x^2*z + la*y^2*x),
  la^2*x^3 + la*y^3 + z^3 + 6*la*x*y*z - t^3 - (z^2*y + la*x^2*z + la*y^2*x)];
polys1 := [&+[Evaluate(MC(f,mon),la3)*mon: mon in mons]: f in polys];
polys1 := Apply(<A1,B1>,polys1);
assert polys1 eq [ x^2*y + la*y^2*z + (la - 1)*z^2*x ,
  la*x*y^2 + la*(la - 1)*y*z^2 + (la - 1)*z*x^2  - (la^2 - la + 1)*t^3 ];

polys := Apply(<A2,B2>,[G1,G2]);
assert polys eq [ x^2*y + x*z^2 + la*y^2*z - t^3,
    x^3 + x^2*z + la*x*y^2 + 6*la*x*y*z + la^2*y^3 + la*y*z^2 
   + la*z^3 + 9*la*t^3 ];
polys1 := [&+[Evaluate(MC(f,mon),la3)*mon: mon in mons]: f in polys];
polys1 := Apply(<A3,B3>,polys1);
assert polys1 eq [ la*(la - 1)*x^2*y + (la - 1)*y^2*z + la*z^2*x,
  (la - 1)*x*y^2 + y*z^2 + la*z*x^2 - (1/(la^2 - la + 1))*t^3 ];

print "Section 5.2 : Minimisation";

R<a,b,p> := PolynomialRing(Rationals(),3);
IsDiv := func<f,subst,g|
  forall{ c : c in Coefficients(Evaluate(f,subst)) | IsDivisibleBy(c,g)}>;

RR<x,y,z> := PolynomialRing(R,3);
F,G := XE7Equations(p*a,p*b:R := RR);
assert IsDiv(G,[x,y,p*z],p^2);
F,G := XE7Equations(p*a,p^2*b:R := RR);
assert IsDiv(F,[p*x,p*y,z],p^3);
assert IsDiv(G,[p*x,y,p*z],p^4);
F,G := XE7Equations(p^2*a,p^2*b:R := RR);
assert IsDiv(F,[x,p*y,z],p^2);
assert IsDiv(G,[p*x,y,p^2*z],p^6);

RR<x,y,z,t> := PolynomialRing(R,4);
F1,F2,G1,G2 := XE9Equations(p*a,p*b:R := RR);
assert IsDiv(F1,[p*x,p*y,p*z,t],p^2) and IsDiv(F2,[p*x,p*y,p*z,t],p^3); 
assert IsDiv(G1,[p*x,y,z,t],p) and IsDiv(G2,[p*x,y,z,t],p);
F1,F2,G1,G2 := XE9Equations(p*a,p^2*b:R := RR);
assert IsDiv(F1,[p*x,p*y,z,t],p^2) and IsDiv(F2,[p*x,p*y,z,t],p^2);
assert IsDiv(G1,[p*x,y,z,t],p) and IsDiv(G2,[p*x,y,z,t],p^2);
F1,F2,G1,G2 := XE9Equations(p^2*a,p^2*b:R := RR);
assert IsDiv(F1,[p^2*x,p^2*y,p*z,t],p^4) and IsDiv(F2,[p^2*x,p^2*y,p*z,t],p^5);
assert IsDiv(G1,[p*x,y,z,t],p^2) and IsDiv(G2,[p*x,y,z,t],p^2);

RR<v,w,x,y,z> := PolynomialRing(R,5);
F,G := XE11Equations(p*a,p*b:R := RR);
assert IsDiv(F,[p*v+w,w,x,y,z],p);
F,G := XE11Equations(p*a,p^2*b:R := RR);
assert IsDiv(F,[p*v,p*w,x,y,z],p^2);
assert IsDiv(G,[p*v,p*w,p*x,y,z],p^2);
F,G := XE11Equations(p^2*a,p^2*b:R := RR);
assert IsDiv(F,[p*v,p*w,x,p*y-x,z],p^3);
assert IsDiv(G,[p^2*v,p*w,p^2*x,p*y,z],p^4);

print "Section 6.1 : Computing the j-invariant";

R<a,b,c,d> := PolynomialRing(Rationals(),4);
RR<r,s> := PolynomialRing(R,2);
F1 := a^2*b + b^2*c + c^2*a;
F2 := a*b^2 + b*c^2 + c*a^2 - d^3;
H1 := Matrix(4,4,[Deriv(Deriv(F1,i),j): i,j in [1..4]]);
H2 := Matrix(4,4,[Deriv(Deriv(F2,i),j): i,j in [1..4]]);
H := Determinant(r*H1+s*H2);
assert H eq 48*(r^3+s^3)*s*(a^3 + b^3 + c^3 - 3*a*b*c)*d;

Ev := Evaluate;
Lambda := func<q,r|Matrix(4,4,[
  [0,Ev(q,[a,b,c,0])*d^2,-Ev(q,[c,a,b,0])*d^2,Ev(r,[a,b,c,0])],
  [-Ev(q,[a,b,c,0])*d^2,0,Ev(q,[b,c,a,0])*d^2,Ev(r,[b,c,a,0])],
  [Ev(q,[c,a,b,0])*d^2,-Ev(q,[b,c,a,0])*d^2,0,Ev(r,[c,a,b,0])],
  [-Ev(r,[a,b,c,0]),-Ev(r,[b,c,a,0]),-Ev(r,[c,a,b,0]),0]])>;
Lambda2 := Lambda(3*(b^2 + 2*c*a),2*(c^3*a + b^3*a) - 3*c^2*b^2 - a^4);
assert Lambda2 eq Matrix(4,4,[&+[MySign(i,j,k,l)*Deriv(F1,l)*Deriv(F2,k):
			       k,l in [1..4]]: i,j in [1..4]]);
D := -(a^3 + b^3 + c^3 - 3*a*b*c)*d;
gamma1 := 18*d^3;  
gamma2 := -6*(a^3 + b^3 + c^3 + 6*a*b*c); 
I := ideal<R|F1,F2>;
MM := Matrix(4,4,[R.i*R.j: i,j in [1..4]]);
assert Eltseq(Lambda2*H1*Lambda2 - gamma1*D*MM) subset I;
assert Eltseq(Lambda2*H2*Lambda2 - gamma2*D*MM) subset I;

invariants := Level9Covariants(R);
_,_,_,c4,c6 := Explode(invariants);
A := d^3;
B := a^3 + b^3 + c^3 + 6*a*b*c;
cc4 := -216*A^3*B + B^4;
cc6 := 5832*A^6 - 540*A^3*B^3 - B^6;
assert [c4-cc4,c6-cc6] subset I;

K<a,b> := FunctionField(Rationals(),2);
R<r,s> := PolynomialRing(K,2);
RR := PolynomialRing(R,4);
F1,F2,G1,G2 := XE9Equations(a,b:R:=RR);

H1 := Matrix(4,4,[Deriv(Deriv(F1,i),j): i,j in [1..4]]);
H2 := Matrix(4,4,[Deriv(Deriv(F2,i),j): i,j in [1..4]]);
H := Determinant(r*H1+s*H2);
f := r^4 + 2*a*r^2*s^2 - 4*b*r*s^3 - (1/3)*a^2*s^4;
D9 := ExactQuotient(H,f);
assert Coefficients(D9) subset K;

la := r;
mu := -3*s;
c4 := -(1/27)*a;
c6 := -(1/54)*b;
DD := la^4 - 6*c4*la^2*mu^2 - 8*c6*la*mu^3 - 3*c4^2*mu^4; 
assert DD eq f;

H1 := Matrix(4,4,[Deriv(Deriv(G1,i),j): i,j in [1..4]]);
H2 := Matrix(4,4,[Deriv(Deriv(G2,i),j): i,j in [1..4]]);
H := Determinant(r*H1+s*H2);
f := a*r^4 - 6*b*r^3*s - 2*a^2*r^2*s^2 + 2*a*b*r*s^3 - ((1/3)*a^3 + 3*b^2)*s^4;
D9 := ExactQuotient(H,f);
assert Coefficients(D9) subset K;

la := r;
mu := -3*s;
c4 := -(1/27)*a;
c6 := -(1/54)*b;
cc4 := c4*la^4 + 4*c6*la^3*mu + 6*c4^2*la^2*mu^2 + 4*c4*c6*la*mu^3 
   + (-3*c4^3 + 4*c6^2)*mu^4;
assert -27*cc4 eq f;

print "Section 6.2 : Modular interpretation of X(n)";

print "Case n = 7";

R<a,b,c> := PolynomialRing(Rationals(),3);
F := a^3*b + b^3*c + c^3*a;
I := ideal<R|F>;
X_7 := Curve(Proj(R),I);
KK := FunctionField(X_7);
la := KK!(-a*c^2/b^3);
tau := KK!(a*c/b^2);
assert tau^7 eq la^4*(la-1);

K<z> := CyclotomicField(7);
e1 := Vector(K,[1,0,0]);
cusps := SetToSequence({X_7(K)!Eltseq(e1*Transpose(g)): g in G7});
assert #cusps eq 24;
X1_7 := ProjectiveSpace(Rationals(),1);
psi := map<X_7->X1_7|[-a*c^2,b^3]>;
KK1<la> := FunctionField(X1_7);
Delta := la*(la - 1)*(la^3 - 8*la^2 + 5*la + 1)^7;
for P in cusps do
  P1 := psi(P);
  assert P1 eq X1_7![1,0] or Evaluate(Delta,P1) eq 0;
end for;

C7 := CC(7:lala:=[la]);
assert Discriminant(C7) eq Delta;
c4C,c6C := Explode(cInvariants(C7));
_,_,c4,c6 := Explode(Level7Invariants(R));

PP := X1_7(KK)![-a*c^2/b^3,1];
xx := [KK|x/c: x in [a,b,c]];
xi7 := KK!(a*c^4/b^5);
assert Evaluate(c4C,PP) eq xi7^4*(&*xx)^2*Evaluate(c4,xx);
assert Evaluate(c6C,PP) eq xi7^6*(&*xx)^3*Evaluate(c6,xx);

print "Case n = 9";

R<a,b,c,d> := PolynomialRing(Rationals(),4);
F1 := a^2*b + b^2*c + c^2*a;
F2 := a*b^2 + b*c^2 + c*a^2 - d^3;
I := ideal<R|F1,F2>;
X_9 := Curve(Proj(R),I);
KK := FunctionField(X_9);

la := KK!(-a*c/b^2);
tau := KK!(a^2*d/(b^2*c));
assert tau^9 eq la*(la-1)^7*(la^2-la+1)^3;

K<z> := CyclotomicField(9);
e1 := Vector(K,[1,0,0,0]);
cusps := SetToSequence({X_9(K)!Eltseq(e1*Transpose(g)): g in G9});
assert #cusps eq 36;
X1_9 := ProjectiveSpace(Rationals(),1);
psi := map<X_9->X1_9|[-a*c,b^2]>;
KK1<la> := FunctionField(X1_9);
Delta := la*(la - 1)*(la^2 - la + 1)^3*(la^3 - 6*la^2 + 3*la + 1)^9;
for P in cusps do
  P1 := psi(P);
  assert P1 eq X1_9![1,0] or Evaluate(Delta,P1) eq 0;
end for;

C9 := CC(9:lala:=[la]);
assert Discriminant(C9) eq Delta;
c4C,c6C := Explode(cInvariants(C9));
invariants := Level9Covariants(R);
_,_,_,c4,c6 := Explode(invariants);

PP := X1_9(KK)![-a*c/b^2,1];
xx := [KK|x/d: x in [a,b,c,d]];
xi9 := KK!(a^2*d^3/(b^4*c));
assert Evaluate(c4C,PP) eq xi9^4*Evaluate(c4,xx);
assert Evaluate(c6C,PP) eq xi9^6*Evaluate(c6,xx);

print "Case n = 11";

R<a,b,c,d,e> := PolynomialRing(Rationals(),5);
F := a^2*b + b^2*c + c^2*d + d^2*e + e^2*a;
H := Matrix(5,5,[Derivative(Derivative(F,i),j): i,j in [1..5]]);
I := ideal<R|Minors(H,4)>;
X_11 := Curve(Proj(R),I);
KK := FunctionField(X_11);

la := KK!(-a*b*d/(c^2*e));
nu := KK!(a*b^3/(c^3*e));
tau := KK!(-a*b/c^2);
assert nu^2 + nu eq la^3 - la^2;
assert tau^11 eq la*nu^2*(la-1)*(la-nu-1)^3;

K<z> := CyclotomicField(11);
e1 := Vector(K,[1,0,0,0,0]);
cusps := SetToSequence({X_11(K)!Eltseq(e1*Transpose(g)): g in G11});
assert #cusps eq 60;
X1_11 := EllipticCurve([0,-1,1,0,0]);
psi := map<X_11->X1_11|[-a*b*c*d,a*b^3,c^3*e]>;
KK1<la,nu> := FunctionField(X1_11);
f := ((-3*la + 2)*nu - la^3 + 5*la^2 - 5*la + 1)/(la-1);
Delta := la*(la-1)*(nu+1)^6*(la*nu + 2*la^2 - 2*la + 1)*f^11;
cusps1 := [c : c in cusps | &*Eltseq(c) ne 0];
assert #cusps1 eq 55;
for P in cusps1 do
  P1 := psi(P);
  assert Evaluate(Delta,P1) eq 0;
end for;

C11 := CC(11:lala:=[la,nu]);
assert Discriminant(C11) eq Delta;
c4C,c6C := Explode(cInvariants(C11));
invariants := Level11Invariants(R);
c4 := invariants[6];
c6tilde := Evaluate(c6tilde,[a,b,c,d,e]);
assert (c6tilde^2 - a*b*c*d*e*(c4^3 - 1728*F^11)) in I;

PP := X1_11(KK)![-a*b*d/(c^2*e),a*b^3/(c^3*e)];
xx := [KK|x/e: x in [a,b,c,d,e]];
xi11 := KK!(a^3*b*e^2/c^6);
assert Evaluate(c4C,PP) eq xi11^4*(&*xx)*Evaluate(c4,xx);
assert Evaluate(c6C,PP) eq xi11^6*(&*xx)*Evaluate(c6tilde,xx);

print "Section 6.3 : An alternative projective embedding";

R<x,y,z> := PolynomialRing(Rationals(),3);
F := x^3*y + y^3*z + z^3*x;
X := Curve(Proj(R),F);
P4<t1,t2,t3,t4> := ProjectiveSpace(Rationals(),3);
psi := map<X->P4|[x^2*y,y^2*z,z^2*x,x*y*z]>;
XX := Image(psi);
M := Matrix(3,4,[t1,0,t4,-t2,t2,-t3,0,t4,t3,t4,-t1,0]);
assert Ideal(Minors(M,3)) eq Ideal(XX);

K0<z7> := CyclotomicField(7);
g7 := 1 + 2*(z7^3 + z7^5 + z7^6);
M2 := (1/g7)*Matrix(K0,3,3,[z7^(i*j)-z7^(-i*j):i,j in [1,2,4]]);
R<x,y,z> := PolynomialRing(K0,3);
F := x^3*y + y^3*z + z^3*x;
f := (7*x*y*z)^M2;
g := x^2*y + y^2*z + z^2*x + x*y*z;
assert IsDivisibleBy(x*y*z*f + g^2,F); 
// => X(7) \cap {g = 0} = 2 T_0 + 2 T_1
 
R<v,w,x,y,z> := PolynomialRing(Rationals(),5);
M5 := Matrix(Rationals(),5,5,[<1,2,1>,<2,3,1>,<3,4,1>,<4,5,1>,<5,1,1>]);
CyclicSum := func<f|&+[f^(M5^i): i in [0..4]]>;
F := CyclicSum(v^2*w);
H := Matrix(5,5,[Derivative(Derivative(F,i),j): i,j in [1..5]]);
I := ideal<R|Minors(H,4)>;
X_11 := Curve(Proj(R),I);
ff := [(v^2*w*x*z)^(M5^i): i in [0..4]];
ff := ff cat [v*w*x*y*z];

P5<t1,t2,t3,t4,t5,t6> := ProjectiveSpace(Rationals(),5);
RR := CoordinateRing(P5);
psi := map<X_11->P5|ff>;
C25 := Image(psi);
assert Degree(C25) eq 25;
h := t6^4 - (t1^2*t2 + t2^2*t3 + t3^2*t4 + t4^2*t5 + t5^2*t1)*t6
  + t1^2*t3*t5 + t2^2*t4*t1 + t3^2*t5*t2 + t4^2*t1*t3 + t5^2*t2*t4;
J := ideal<RR|[Deriv(h,i): i in [1..6]]>;
assert J subset Ideal(C25);
assert Ideal(C25)^2 subset J;

V := VectorSpace(Rationals(),5);
PP := [X_11!Eltseq(V.i): i in [1..5]];
KK := FunctionField(X_11);
// T0 := &+[Divisor(P): P in PP]; 
// assert Dimension(RiemannRochSpace(5*T0)) eq 6;

// print "Checking P1 + 6 P_3 + 3 P_4 + 10 P_5 \sim H";
// nn := [1,0,6,3,10];
// for i in [1..5] do
//  time assert Valuation(KK!(R.2/R.i),PP[i]) eq nn[i];
// end for;

K0<z11> := CyclotomicField(11);
g11 := 1 + 2*&+[z11^i: i in [1,3,9,5,4]];
M2 := (1/g11)*Matrix(K0,5,5,[z11^(i*j)-z11^(-i*j):i,j in [1,3,9,5,4]]);
M5 := Matrix(K0,5,5,[<1,2,1>,<2,3,1>,<3,4,1>,<4,5,1>,<5,1,1>]);
R<v,w,x,y,z> := PolynomialRing(K0,5);
F := v^2*w + w^2*x + x^2*y + y^2*z + z^2*v;
H := Matrix(5,5,[Derivative(Derivative(F,i),j): i,j in [1..5]]);
I := ideal<R|Minors(H,4)>;
f := (v*w*x*y*z)^M2;
g := &+[(v^2*w*x*z)^(M5^i): i in [0..4]] - v*w*x*y*z;
assert 121*f*(v*w*x*y*z)^3 - g^4 in I; 
// => X(11) \cap {g = 0} = 15 T_0 + 5 T_1 

print "Section 6.4 : Formulae in the case n = 7";

R<a,b,c> := PolynomialRing(Rationals(),3);
RR<x,y,z> := PolynomialRing(R,3);

F := x^3*y + y^3*z + z^3*x;

s1 := (a^2*c^3 - 2*a*b^3*c)*x^2*y + (a^3*b^2 - 2*a*b*c^3)*y^2*z 
    + (b^3*c^2 - 2*a^3*b*c)*x*z^2 + (a^3*c^2 + a^2*b^3 + b^2*c^3)*x*y*z;
s2 := a^2*b*x^2*y + b^2*c*y^2*z + c^2*a*z^2*x + 2*a*b*c*x*y*z;

mons3 := MD(R,3);
mons6 := MD(R,6);
I := ideal<R|a^3*b + b^3*c + c^3*a>;
ff := [NormalForm(mon*a*b*c,I): mon in mons3];
mat1 := Matrix(#ff,#mons6,[MC(f,mon):mon in mons6,f in ff]);
mons := Monomials(s1^2);
gg := [NormalForm(R!mon,I) : mon in mons];
mat2 := Matrix(#gg,#mons6,[MC(g,mon):mon in mons6,g in gg]);
mat := Solution(mat1,mat2);
mons := [Evaluate(mon,[x,y,z]): mon in mons];
mons3 := [Evaluate(mon,[x,y,z]): mon in mons3];
r1 := &+[&+[MC(s1^2,mons[i])*mat[i,j]: i in [1..10]]*mons3[j]: j in [1..10]];
r2 := &+[&+[MC(s2^2,mons[i])*mat[i,j]: i in [1..10]]*mons3[j]: j in [1..10]];
assert IsDivisibleBy(r1*x*y*z - s1^2,F);
assert IsDivisibleBy(r2*x*y*z - s2^2,F);

invariants := Level7Invariants(R);
F,H,c4,c6 := Explode(invariants);
assert c4^3 - c6^2 eq 1728*H^7 - 2048*F^9*H + 256*F^7*c4 + 22016*F^6*H^3 
  - 1088*F^4*H^2*c4 - 60032*F^3*H^5 - 88*F^2*H*c4^2 - 1008*F*H^4*c4;
Delta := H^7;
Nabla := func<f|Vector(3,[Deriv(f,i): i in [1..3]])>;
xx,ee,ff := Explode(Level7Covariants(R));

FF,GG := XE7EquationsRaw(c4,c6:R:=RR);
d1 := -6*(3*x^2 + c4*x*z - 3*c4*y^2 + c6*y*z)*z; 
d2 := 2*Delta*(4*x^3 + c4*x^2*z - 12*c4*x*y^2 - 2*c6*x*y*z 
	       + 8*c6*y^3 + c4^2*y^2*z + 200*Delta*z^3);
I := ideal<R|a^3*b + b^3*c + c^3*a>;

covs := [xx,Cross(Nabla(F),Nabla(H)),H*ee];
subst1 := [&+[covs[i][j]*RR.i: i in [1..3]]: j in [1..3]];
diffr := H^3*FF - Evaluate(F,subst1);
assert Coefficients(diffr) subset I;
diffr := 2*a*b*c*H^4*d1 - Evaluate(r1,subst1);
assert Coefficients(diffr) subset I;

cons := [Nabla(F),Cross(xx,ee),H^2*Nabla(H)];
subst2 := [&+[cons[i][j]*RR.i: i in [1..3]]: j in [1..3]];
diffr := H^2*GG - Evaluate(F,subst2);
assert Coefficients(diffr) subset I;
diffr := a*b*c*d2 - 2*H^5*Evaluate(r2,subst2);
assert Coefficients(diffr) subset I;

K<a,b> := FunctionField(Rationals(),2);
R0 := PolynomialRing(Integers(),2);
R<x,y,z> := PolynomialRing(K,3);

d_11 := -2*(a*x^2 + 3*b*x*z + 3*y^2 + 2*a*y*z)*z;
d_12 := 2*(a*x^2 + 3*b*x*z + 3*y^2 + 2*a*y*z)*x;
d_13 := 4*(3*b*x^2 - 2*a*x*y - 2*a^2*x*z - 3*b*y*z - 2*a*b*z^2)*z;
d_14 := 4*(a^2*x^2 + 3*b*x*y + 4*a*b*x*z + a*y^2 + 3*b^2*z^2)*z;

dd_11 := -7*a*x^2*y + 6*x^2*z + 3*a^2*y^3 - 8*a*y^2*z + 3*y*z^2;
dd_12 := 2*a*x^3 + 12*b*x^2*y - 2*a*x*y*z - 3*a*b*y^3 + 6*b*y^2*z;
dd_13 := 2*a^2*x*y^2 - 10*a*x*y*z + 6*x*z^2 + 5*a*b*y^3 - 12*b*y^2*z;
dd_14 := 2*a^2*x^2*y - 3*a*x^2*z + 5*a*b*x*y^2 - 12*b*x*y*z
             - 3*a^2*y^2*z + 8*a*y*z^2 - 3*z^3;

F,G := XE7Equations(a,b:R:=R);
c4 := -(1/27)*a;
c6 := -(1/54)*b;
Delta := (c4^3 - c6^2)/1728;
d1 := -6*(3*x^2 + c4*x*z - 3*c4*y^2 + c6*y*z)*z; 
d2 := 2*Delta*(4*x^3 + c4*x^2*z - 12*c4*x*y^2 - 2*c6*x*y*z 
	       + 8*c6*y^3 + c4^2*y^2*z + 200*Delta*z^3);
subst := [6*c4*z-(1/3)*y,x,-18*z];
mat := Matrix(3,3,[MC(f,R.j): j in [1..3],f in subst]);
d1 := Determinant(mat)*Evaluate(d1,subst);
subst := [9*c4*y+z,3*x,108*y];
mat := Matrix(3,3,[MC(f,R.j): j in [1..3],f in subst]);
d2 := Determinant(mat)*Evaluate(d2,subst);

assert d1 eq 6^2*d_11;

I := IdealWithFixedBasis([d_11,F]);
d := [d_11,d_12,d_13,d_14];
coeffs := Matrix(4,10,[MC(f,mon): mon in MD(R,3),f in d]);
assert K!GCD([R0!Numerator(m): m in Minors(coeffs,4)]) eq 2^6*(4*a^3+27*b^2);
dd := Matrix(4,4,[Coordinates(I,d[i]*d[j])[1]:i,j in [1..4]]);
assert Minors(dd,2) subset ideal<R|F>;

assert [dd[i,i]: i in [1..4]] eq
  [-2*a*x^2*z - 6*b*x*z^2 - 6*y^2*z - 4*a*y*z^2,
  8*b*x^3 - 4*a*x^2*y - 6*a^2*x^2*z - 12*b*x*y*z - 10*a*b*x*z^2 + 4*y^3 + 
  6*a*y^2*z + 4*a^2*y*z^2 - 8*b^2*z^3,
  18*b*x^3 - 12*a*x^2*y - 8*a^2*x^2*z + 12*b*x*y*z + 10*a*b*x*z^2 - 
  8*a*y^2*z - 12*a^2*y*z^2 + (-8*a^3 - 24*b^2)*z^3,
  -6*a*b*x^3 + 4*a^2*x^2*y + (-8*a^3 - 42*b^2)*x^2*z - 20*a*b*x*y*z - 
  22*a^2*b*x*z^2 + (4*a^3 - 12*b^2)*y*z^2 - 26*a*b^2*z^3];

XE7 := Curve(Proj(R),F);
KK<u,v> := FunctionField(XE7);
ff := KK!(d1/y^3);
vv,ww := Support(Divisor(ff));
D1 := &+[Integers()!((1/2)*ww[i])*vv[i]: i in [1..#vv] | ww[i] gt 0];
ff := KK!(z/y);
vv,ww := Support(Divisor(ff));
H := &+[ww[i]*vv[i]: i in [1..#vv] | ww[i] gt 0];
assert Degree(D1) eq 6;
assert Degree(H) eq 4;
V,Vmap := RiemannRochSpace(3*H - D1);
assert Dimension(V) eq 4;
M := Matrix(4,4,[Eltseq((KK!(dd[1,i]/z^3)) @@ Vmap): i in [1..4]]);
assert Determinant(M) ne 0;

S<t1,t2,t3,t4> := PolynomialRing(K,4);
M := Matrix(3,4,[ 0, t3 , t4 , -2*a*t1 - t4, t1 , -2*a*t1 - t4, 
  2*b*t1 + a*t2 + a*t3 , 2*a*t2 + a*t3, t2, -2*b*t1 - a*t3, 
  -a^2*t1 + b*t3 - a*t4 , 2*b*t2 - b*t3 - a*t4 ]);
psi := map<XE7->Proj(S)|[d_11,d_12,d_13,d_14]>;
assert ideal<S|Minors(M,3)> eq Ideal(Image(psi));

I := IdealWithFixedBasis([dd_11,G]);
d := [dd_11,dd_12,dd_13,dd_14];
coeffs := Matrix(4,10,[MC(f,mon): mon in MD(R,3),f in d]);
assert K!GCD([R0!Numerator(m): m in Minors(coeffs,4)]) eq 2;
dd := Matrix(4,4,[Coordinates(I,d[i]*d[j])[1]:i,j in [1..4]]);
assert Minors(dd,2) subset ideal<R|G>;
assert [dd[i,i]: i in [1..4]] eq
  [-7*a*x^2*y + 6*x^2*z + 3*a^2*y^3 - 8*a*y^2*z + 3*y*z^2,
  -8*b*x^3 + 3*a^2*x^2*y + 4*a*x^2*z - 2*a*b*x*y^2 + 12*b*x*y*z 
    + (-3*a^3 - 16*b^2)*y^3 + 2*a^2*y^2*z - 3*a*y*z^2 + 2*z^3,
  -a^2*x^2*y + 2*a*b*x*y^2 - 12*b*x*y*z + (-7*a^3 - 36*b^2)*y^3 
    + 6*a^2*y^2*z - 7*a*y*z^2 + 6*z^3,
  -a^3*x^2*y + 2*a^2*x^2*z - 2*a^2*b*x*y^2 + 2*a*b*x*y*z + 6*b*x*z^2 
    + (-3*a^4 - 19*a*b^2)*y^3 + (8*a^3 + 42*b^2)*y^2*z - 3*a^2*y*z^2];

E := EllipticCurve([a,b]);
Delta := Discriminant(E);
assert (2^5*3^4)^2*d2 eq Delta*(dd[2,2] + 2*dd[2,3] + dd[3,3]); 

XEm7 := Curve(Proj(R),G);
KK<u,v> := FunctionField(XEm7);
ff := KK!(dd_11/y^3);
vv,ww := Support(Divisor(ff));
D1 := &+[Integers()!((1/2)*ww[i])*vv[i]: i in [1..#vv] | ww[i] gt 0];
ff := KK!(z/y);
vv,ww := Support(Divisor(ff));
H := &+[ww[i]*vv[i]: i in [1..#vv] | ww[i] gt 0];
assert Degree(D1) eq 6;
assert Degree(H) eq 4;
V,Vmap := RiemannRochSpace(3*H - D1);
assert Dimension(V) eq 4;
M := Matrix(4,4,[Eltseq((KK!(dd[1,i]/z^3)) @@ Vmap): i in [1..4]]);
assert Determinant(M) ne 0;

print "Section 6.5 : Formulae in the case n = 11";

function MyEvaluate(poly,subst,I)
  RR := Parent(poly);
  indices := [[0,0,0,0,0]];
  prods := [RR!1];
  ans := 0;
  mmm := MD(RR,5);
  for mon in mmm do	 
    coeff := NormalForm(MC(poly,mon),I);
    printf ".";
    ii := [Degree(mon,j): j in [1..5]];
    rr := [jj : jj in indices | forall{i : i in [1..5] | ii[i] ge jj[i]}];
    _,t := Maximum([&+r: r in rr]);
    jj := rr[t];
    aa := Position(indices,jj);
    pp := prods[aa];
    for i in [1..5] do
      ee := ii[i] - jj[i];
      assert ee ge 0;
      for k in [1..ee] do
        pp := pp*subst[i];
        pp := &+[NormalForm(Coefficients(pp)[j],I)*Monomials(pp)[j]
                                             : j in [1..#Terms(pp)]];
        jj[i] +:= 1;
        Append(~indices,jj);
        Append(~prods,pp);
      end for;
    end for;
    ans +:= coeff*pp;
  end for;
  print "";
  return ans;
end function;

function s1(aa,xx)
  a,b,c,d,e := Explode(aa);
  v,w,x,y,z := Explode(xx);
  return (a^3*b*c^3 + b^4*c*d^2 - a*b^2*c^2*d*e - 2*b*c^2*d*e^3)*v^2*w*x*z
    + (b^3*c*d^3 + c^4*d*e^2 - a*b*c^2*d^2*e - 2*a^3*c*d^2*e)*v*w^2*x*y
    + (c^3*d*e^3 + a^2*d^4*e - a*b*c*d^2*e^2 - 2*a*b^3*d*e^2)*w*x^2*y*z
    + (a^3*d^3*e + a*b^2*e^4 - a^2*b*c*d*e^2 - 2*a^2*b*c^3*e)*v*x*y^2*z
    + (a*b^3*e^3 + a^4*b*c^2 - a^2*b^2*c*d*e - 2*a*b^2*c*d^3)*v*w*y*z^2
    + 2*(a^2*b^2*c^2*e + a^2*b^2*d*e^2 + a^2*c*d^2*e^2 + a*b^2*c^2*d^2
	+ b*c^2*d^2*e^2)*v*w*x*y*z;
end function;

function s2(aa,xx)
  a,b,c,d,e := Explode(aa);
  v,w,x,y,z := Explode(xx);
  return a^2*b*c*e*v^2*w*x*z + a*b^2*c*d*v*w^2*x*y + b*c^2*d*e*w*x^2*y*z 
    + a*c*d^2*e*v*x*y^2*z + a*b*d*e^2*v*w*y*z^2 + 2*a*b*c*d*e*v*w*x*y*z;
end function;

load "quintics.m";

R<a,b,c,d,e> := PolynomialRing(Rationals(),5);
F,H,I7,I8,I9,c4 := Explode(Level11Invariants(R));
Hmat := Matrix(5,5,[Derivative(Derivative(F,i),j): i,j in [1..5]]);
I := ideal<R|Minors(Hmat,4)>;

K := FieldOfFractions(R); 
RR<v,w,x,y,z> := PolynomialRing(K,5);
aa := [a,b,c,d,e];
xx := [v,w,x,y,z];
II := ideal<RR|[Evaluate(f,xx): f in Basis(I)]>;
nf := NormalForm((&*aa)^3*(&*xx)^3*r1tilde(aa,xx) - s1(aa,xx)^4,II);
assert forall{ f : f in Coefficients(nf) | R!f in I};
nf := NormalForm((&*aa)^3*(&*xx)^3*r2tilde(aa,xx) - s2(aa,xx)^4,II);
assert forall{ f : f in Coefficients(nf) | R!f in I};

RR<v,w,x,y,z> := PolynomialRing(R,5);
Nabla := func<f|Vector(5,[Deriv(f,i): i in [1..5]])>;
covs := Level11Covariants(R);
cons := [Nabla(f): f in [F,I7,I8,I9,c4]];

subst1 := [&+[covs[i][j]*RR.i: i in [1..5]]: j in [1..5]];
subst2 := [&+[cons[i][j]*RR.i: i in [1..5]]: j in [1..5]];

print "Twisting the quintic forms";

xx := [v,w,x,y,z];
xx0 := [v,w,x,0,0]; 

// For a complete verification of the formulae the following line
// should be uncommented.
// xx0 := [v,w,x,y,z]; // Warning : Takes a few hours.

r1t := r1tilde(aa,xx);
subst1 := [Evaluate(f,xx0): f in subst1]; 
time d1new := MyEvaluate(r1t,subst1,I);
// time d1old := Evaluate(r1t,subst1);
// time assert Coefficients(d1new-d1old) subset I;

r2t := r2tilde(aa,xx);
subst2 := [Evaluate(f,xx0): f in subst2]; 
time d2new := MyEvaluate(r2t,subst2,I);
// time d2old := Evaluate(r2t,subst2);
// time assert Coefficients(d2new-d2old) subset I;

inv := [F,c4];
S<F,c4> := PolynomialRing(Rationals(),2);
SS<v,w,x,y,z> := PolynomialRing(S,5);
xx := [v,w,x,y,z];

d1t := d1tilde(F,c4,xx);
d2t := d2tilde(F,c4,xx);

time for mon in MD(RR,5) do
  m := Evaluate(mon,xx0);
  if m ne 0 then 
    printf ".";
    assert MC(d1new,m) - Evaluate(MC(d1t,SS!m),inv) in I;
    assert MC(d2new,m) - Evaluate(MC(d2t,SS!m),inv) in I;
  end if;
end for;
print "";

subst := [F*v,w,F^7*x,F^2*y,F^4*z];
assert d1(c4,F^11,xx) eq Evaluate(d1t,subst);

subst := [F^2*v,F^8*w,F^4*x,y,F^3*z];
assert F^4*d2(c4,F^11,xx) eq Evaluate(d2t,subst);

print "2-isogenous curves";

K<a,b> := FunctionField(Rationals(),2);
R<v,w,x,y,z> := PolynomialRing(K,5);
_,G := XE11Equations(a,b:R:=R);
H := Matrix(5,5,[Deriv(Deriv(G,i),j):i,j in [1..5]]);
X := Scheme(Proj(R),Minors(H,4));
P<T> := PolynomialRing(K);
L<t> := quo<P|T^3 + a*T + b>; 
P := X(L)![0,-t^2 - a,1/2*(t^2 + a),t,1];
c4 := -a/27;
c6 := -b/54;
Delta := (c4^3-c6^2)/1728;
_,GG := XE11EquationsRaw(c4,Delta:R:=R);
v11 := 44*(2*c4*v - 6*c6*w + 33*c6*x + 135*c4^2*y + 810*c4*c6*z);
x11 := 60*(5*v + 729*c4*y + 2187*c6*z);
y11 := 11*(c4*v - 3*c6*w - 6*c6*x);
z11 := 60*(v + 27*c4*y + 81*c6*z);
subst2 := [v11,-427680*y,x11,-y11,-z11];
HH := Matrix(5,5,[Deriv(Deriv(GG,i),j):i,j in [1..5]]);
X := Scheme(Proj(R),Minors(HH,4));
P := X(L)![Evaluate(f,Eltseq(P)): f in subst2];
assert d2(c4,Delta,Eltseq(P)) eq 0;

K<T> := FunctionField(Rationals());
RR<v,w,x,y,z> := PolynomialRing(K,5);
f5 := (T^2 + 10*T + 5)^3;
f7 := (T^2 + 13*T + 49)*(T^2 + 5*T + 1)^3;
f13 := (T^2 + 5*T + 13)*(T^4 + 7*T^3 + 20*T^2 + 19*T + 1)^3;

print "5-isogenous curves";

a := -3*T^4 - 96*T^3 - 1050*T^2 - 4080*T - 1875;
b := -2*T^6 - 96*T^5 - 1818*T^4 - 16784*T^3 - 73782*T^2 - 114000*T + 31250;
E := EllipticCurve([a,b]);
assert jInvariant(E) eq f5/T;
F := XE11Equations(a,b:R := RR);
H := Matrix(5,5,[Deriv(Deriv(F,i),j):i,j in [1..5]]);
X := Scheme(Proj(RR),Minors(H,4));
pol := (T^2 + 22*T + 125)*(T^3 + 15*T^2 + 7*T + 25);
P := X![pol,-pol,-(T - 1)*(T^2 + 22*T + 125),0,T+5];
c4 := -a/27;
c6 := -b/54;
Delta := (c4^3-c6^2)/1728;
v1 := c6*v + 2*c6*w - 6*c4^2*x + 3*c4^2*y - 9*c4*c6*z;
w1 := c6*v + 6*c4^2*x + 3*c4^2*y + 9*c4*c6*z;
subst1 := [-v1,w1,-864*x,-36*c4*x - 108*c6*z,72*y];
P := [Evaluate(f,Eltseq(P)): f in subst1];
assert d1(c4,Delta,[1,0,0,0,0]) eq 0;
assert d1(c4,Delta,P) eq 0;

print "7-isogenous curves";

a := -3*T^4 - 54*T^3 - 345*T^2 - 774*T - 147;
b := -2*T^6 - 54*T^5 - 588*T^4 - 3150*T^3 - 7980*T^2 - 6678*T + 686;
E := EllipticCurve([a,b]);
assert jInvariant(E) eq f7/T; 
_,G := XE11Equations(a,b:R := RR);
H := Matrix(5,5,[Deriv(Deriv(G,i),j):i,j in [1..5]]);
X := Scheme(Proj(RR),Minors(H,4));
P := X![ -36*(T^2 + 9*T - 7)*(T^2 + 13*T + 49),
     2*(T^3 + 10*T^2 + 26*T - 1)*(T^2 + 13*T + 49),
      -(T + 1)*(T^2 + 9*T + 17)*(T^2 + 13*T + 49),
                    -(T + 1)*(T^2 + 13*T + 49),T+5];
c4 := -a/27;
c6 := -b/54;
Delta := (c4^3-c6^2)/1728;
v11 := 44*(2*c4*v - 6*c6*w + 33*c6*x + 135*c4^2*y + 810*c4*c6*z);
x11 := 60*(5*v + 729*c4*y + 2187*c6*z);
y11 := 11*(c4*v - 3*c6*w - 6*c6*x);
z11 := 60*(v + 27*c4*y + 81*c6*z);
subst2 := [v11,-427680*y,x11,-y11,-z11];
P := [Evaluate(f,Eltseq(P)): f in subst2];
assert d2(c4,Delta,P) eq 0;

print "13-isogenous curves";

a := -3*T^8 - 54*T^7 - 459*T^6 - 2322*T^5 - 7500*T^4 - 15354*T^3 
  - 18459*T^2 - 10062*T - 507;
b := -2*T^12 - 54*T^11 - 702*T^10 - 5724*T^9 - 32286*T^8 - 131454*T^7 
  - 391896*T^6 - 849330*T^5 - 1297914*T^4 - 1307700*T^3 - 750906*T^2 
  - 161226*T + 4394;
E := EllipticCurve([a,b]);
assert jInvariant(E) eq f13/T;
_,G := XE11Equations(a,b:R := RR);
H := Matrix(5,5,[Deriv(Deriv(G,i),j):i,j in [1..5]]);
X := Scheme(Proj(RR),Minors(H,4));
P := X![36*(T^2 + 5*T + 13)*(T^2 + 6*T + 13)^2*(T^3 + 4*T^2 + 8*T - 1),
  2*(T^2 + 5*T + 13)*(T^2 + 6*T + 13)*(T^6 + 12*T^5 + 65*T^4 
                                   + 189*T^3 + 296*T^2 + 201*T - 8),
  -(T + 2)*(T^2 + 5*T + 13)*(T^2 + 6*T + 13)^2*(T^3 + 4*T^2 + 8*T - 1),
  -(T + 1)*(T + 2)*(T^2 + 5*T + 13)*(T^2 + 6*T + 13),T^2 + 5*T + 10];
c4 := -a/27;
c6 := -b/54;
Delta := (c4^3-c6^2)/1728;
v11 := 44*(2*c4*v - 6*c6*w + 33*c6*x + 135*c4^2*y + 810*c4*c6*z);
x11 := 60*(5*v + 729*c4*y + 2187*c6*z);
y11 := 11*(c4*v - 3*c6*w - 6*c6*x);
z11 := 60*(v + 27*c4*y + 81*c6*z);
subst2 := [v11,-427680*y,x11,-y11,-z11];
P := [Evaluate(f,Eltseq(P)): f in subst2];
assert d2(c4,Delta,P) eq 0;

print "Section 7 : Examples";

CD := CremonaDatabase();
CR := func<E|Conductor(E) lt LargestConductor(CD) select 
    CremonaReference(E) else IntegerToString(Conductor(E)) cat "*">;

print "Example 7.1";

E := EllipticCurve("162c1");
a := 3645;
b := -13122;
assert IsIsomorphic(EllipticCurve([a,b]),E);
R<x,y,z> := PolynomialRing(Rationals(),3);
F,G := XE7Equations(a,b:R := R);
assert F eq 3645*x^4 - 91854*x^3*z + 3*x^2*y^2 - 39858075*x^2*z^2
  + 78732*x*y*z^2 + 239148450*x*z^3 + 2*y^3*z + 10935*y^2*z^2 
  + 26572050*y*z^3 - 688747536*z^4;
inv :=  -4*(4*a^3+27*b^2);
assert inv eq -2^11*3^18;
mu := 2^10*3^14;
subst1 := [ 36*y - 9*z, 1944*x - 972*y - 1215*z, z ];
FF := (1/mu)*Evaluate(F,subst1);
M := Matrix(Rationals(),3,3,[MC(f,R.i): i in [1..3],f in subst1]);
assert FF eq 3*x^3*z + 3*x^2*y^2 - 6*x^2*y*z + 3*x^2*z^2 - 3*x*y^3 
  + 3*x*z^3 + 2*y^4 - y^3*z - 9*y^2*z^2 + 4*y*z^3 - 5*z^4;
inv := (1/mu)^3*Determinant(M)^4*inv;
assert inv eq -2*3^4;

assert G eq -13286025*x^4 - 95659380*x^3*y + 157464*x^3*z 
  - 296764094574*x^2*y^2 + 21870*x^2*z^2 - 348678440100*x*y^3 
  + 573956280*x*y^2*z - 236196*x*y*z^2 + 541480183553295*y^4 
  - 394652338128*y^3*z + 79716150*y^2*z^2 - 29160*y*z^3 + 3*z^4;
inv := 16*(4*a^3+27*b^2)^2;
assert inv eq 2^22*3^36;
mu := 2^12*3^20;
subst2 := [18*x + 18*y + 9*z,z,-486*x + 1458*y + 1944*z];
GG := (1/mu)*Evaluate(G,subst2);
M := Matrix(Rationals(),3,3,[MC(f,R.i): i in [1..3],f in subst2]);
assert GG eq -x^3*y - x^3*z - 6*x^2*z^2 + 6*x*y^2*z - 6*x*y*z^2 
  + 6*x*z^3 + 2*y^4 + 2*y^3*z - 6*y^2*z^2 - 38*y*z^3 - 8*z^4;
inv := (1/mu)^3*Determinant(M)^4*inv;
assert inv eq 2^2*3^4;

cubics1 := 
  [-2*a*x^2*z - 6*b*x*z^2 - 6*y^2*z - 4*a*y*z^2,
  8*b*x^3 - 4*a*x^2*y - 6*a^2*x^2*z - 12*b*x*y*z - 10*a*b*x*z^2 + 4*y^3 + 
  6*a*y^2*z + 4*a^2*y*z^2 - 8*b^2*z^3,
  18*b*x^3 - 12*a*x^2*y - 8*a^2*x^2*z + 12*b*x*y*z + 10*a*b*x*z^2 - 
  8*a*y^2*z - 12*a^2*y*z^2 + (-8*a^3 - 24*b^2)*z^3,
  -6*a*b*x^3 + 4*a^2*x^2*y + (-8*a^3 - 42*b^2)*x^2*z - 20*a*b*x*y*z - 
  22*a^2*b*x*z^2 + (4*a^3 - 12*b^2)*y*z^2 - 26*a*b^2*z^3];

cubics2 := 
  [-7*a*x^2*y + 6*x^2*z + 3*a^2*y^3 - 8*a*y^2*z + 3*y*z^2,
  -8*b*x^3 + 3*a^2*x^2*y + 4*a*x^2*z - 2*a*b*x*y^2 + 12*b*x*y*z 
    + (-3*a^3 - 16*b^2)*y^3 + 2*a^2*y^2*z - 3*a*y*z^2 + 2*z^3,
  -a^2*x^2*y + 2*a*b*x*y^2 - 12*b*x*y*z + (-7*a^3 - 36*b^2)*y^3 
    + 6*a^2*y^2*z - 7*a*y*z^2 + 6*z^3,
  -a^3*x^2*y + 2*a^2*x^2*z - 2*a^2*b*x*y^2 + 2*a*b*x*y*z + 6*b*x*z^2 
    + (-3*a^4 - 19*a*b^2)*y^3 + (8*a^3 + 42*b^2)*y^2*z - 3*a^2*y*z^2];

function Covariants(F)
  Deriv := Derivative;
  Hmat := Matrix(3,3,[Deriv(Deriv(F,i),j): i,j in [1..3]]);
  Psi_6 := (-1/54)*Determinant(Hmat);
  vv := Matrix(3,1,[Deriv(Psi_6,i): i in [1..3]]);
  HHmat := HorizontalJoin(Hmat,vv);
  HHmat := VerticalJoin(HHmat,Matrix(1,4,Eltseq(vv) cat [0]));
  Psi_14 := (1/9)*Determinant(HHmat);
  Jmat := Matrix(3,3,[Deriv(f,i): i in [1..3],f in [F,Psi_6,Psi_14]]);
  Psi_21 := (1/14)*Determinant(Jmat);
  return Psi_6,Psi_14,Psi_21;
end function;

H,cc4,cc6 := Covariants(F);
pts := [[1,0,0],[3,-2,-1],[1,0,0],[1,1,-1],[4,-1,1]];
subst := subst1;
cubics := cubics1;
ans := ["162c1","293706*","162c2","17334f1","624186*"];

for i in [1..#pts] do
  P := pts[i];
  assert Evaluate(FF,P) eq 0;
  P1 := [Evaluate(f,P): f in subst];
  assert Evaluate(F,P1) eq 0;
  c4 := Evaluate(cc4,P1);
  c6 := Evaluate(cc6,P1);
  assert exists(d){Evaluate(dd,P1): dd  in cubics | Evaluate(dd,P1) ne 0};
  E1 := MinimalModel(EllipticCurve([-27*d^2*c4,-54*d^3*c6]));
  print P,CR(E1),aInvariants(E1);
  if i eq 2 then 
    F := G; FF := GG;
    H,cc4,cc6 := Covariants(F);
    subst := subst2;
    cubics := [Discriminant(E)*f : f in cubics2];
  end if;
  assert CR(E1) eq ans[i];
end for;

print "Example 7.2";

QT<T> := FunctionField(Rationals());
j := 27*T^3*(5*T - 56)/(T - 1);
a := -27*j/(4*(j - 1728));
R := PolynomialRing(Rationals(),3);
RR<x,y,z> := PolynomialRing(QT,3);
_,G := XE7Equations(a,a:R := RR);
P := [0,-4*(T^2-12*T+8)*(5*T^2+4*T+8),9*T^2*(T+4)*(5*T-56)];
assert Evaluate(G,P) eq 0;
Tlist := [-16,8,2,16/5];
dlist := [-38,-10,-2,-42];
cubic := -7*a*x^2*y + 6*x^2*z + 3*a^2*y^3 - 8*a*y^2*z + 3*y*z^2;
cubic := Evaluate(cubic,P);
for i in [1..4] do
  aa := Evaluate(a,Tlist[i]);
  E1 := EllipticCurve([aa,aa]);
  E1 := MinimalModel(QuadraticTwist(E1,dlist[i]));
  GG := &+[Evaluate(Coefficients(G)[j],Tlist[i])*(R!Monomials(G)[j])
    : j in [1..#Terms(G)]];
  H,cc4,cc6 := Covariants(GG);
  c4 := Evaluate(Evaluate(cc4,P),Tlist[i]);
  c6 := Evaluate(Evaluate(cc6,P),Tlist[i]);
  dd := Discriminant(E1)*Evaluate(cubic,Tlist[i]);
  E2 := EllipticCurve([-27*dd^2*c4,-54*dd^3*c6]);
  E2 := MinimalModel(QuadraticTwist(E2,dlist[i]));
  print Tlist[i],dlist[i],CR(E1),CR(E2);
end for;

print "Example 7.3";

E := EllipticCurve("47775z1");
a := -41489280;
b := 102867483600;
assert IsIsomorphic(EllipticCurve([a,b]),E);
R<x,y,z,t> := PolynomialRing(Rationals(),4);
F1,F2,G1,G2 := XE9Equations(a,b:R := R);
inv :=  -2^10*(4*a^3+27*b^2)^4;
assert inv eq -2^42*3^60*5^12*7^32*13^4;
M := Matrix(Rationals(),4,4,[2520473760, 937149484320, 
  -1998984627360, -152410870080, 0, 79644600, -185343480, 
  -3827880, 0, -22932, 47040, 6468, 0, -6, 13, 1]);
subst := [&+[M[i,j]*R.j: j in [1..4]]:i in [1..4]];
rr := [Rationals()|-2^10*3^18*5^2*7^8,-2^12*3^19*5^2*7^10];
N := Matrix(Rationals(),2,2,[13,-2,28,-5])*DiagonalMatrix(rr)^(-1);
FF := [&+[N[i,j]*[F1,F2][j]^M: j in [1..2]]: i in [1..2]];

assert FF eq [ -x^2*z + x^2*t + 4*x*y*z + 2*x*y*t - 3*x*z^2 + 2*x*z*t 
  - 3*x*t^2  + 6*y^3 + 14*y^2*z + y^2*t + 6*y*z^2 - 4*y*z*t + 9*y*t^2 
  - 6*z^3 + 27*z^2*t - 13*z*t^2 - t^3,
  -3*x^2*y + 4*x^2*z + 3*x^2*t + 3*x*y^2 + 20*x*y*z - 12*x*y*t - 3*x*z^2  
  - 32*x*z*t + 25*x*t^2  + 21*y^3 + 16*y^2*z - 24*y^2*t - 12*y*z^2 
  + 100*y*z*t + 34*y*t^2 + 39*z^3 - 21*z^2*t - 56*z*t^2 - 11*t^3];
inv := Determinant(N)^6*Determinant(M)^9*inv;
assert inv eq -3^3*5^6*7^5*13^4;

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

pts := [[1,0,0,0],[4,-1,-1,0],[1,2,-1,0]];
c4 := -a/27;
c6 := -b/54;

for i in [1..3] do
  P := pts[i];
  assert Evaluate(FF[1],P) eq 0 and Evaluate(FF[2],P) eq 0;
  P1 := [Evaluate(f,P): f in subst];
  assert Evaluate(F1,P1) eq 0 and Evaluate(F2,P1) eq 0;
  gamma1,gamma2 := Explode(TangentLineConstruction([F1,F2],P1));
  la := gamma2;
  mu := 3*gamma1;
  _,cc4,cc6 := HessePolynomials(3,1,[c4,c6]:Variables := [la,mu]);
  E1 := MinimalModel(EllipticCurve([-27*cc4,-54*cc6]));
  print CR(E1),aInvariants(E1);
end for;

X := Curve(Proj(R),[G1,G2]);
// assert not IsLocallySoluble(X,7);

print "Example 7.4";

E := EllipticCurve("201c1");
a := -1029699;
b := 402173694;
assert IsIsomorphic(EllipticCurve([a,b]),E);
R<x,y,z,t> := PolynomialRing(Rationals(),4);
F1,F2,G1,G2 := XE9Equations(a,b:R := R);
inv :=  2^8*(4*a^3+27*b^2)^5;
assert inv eq 2^48*3^85*67^5;
M := Matrix(Rationals(),4,4,[ -26471709, -23136696, 20106774, 
  -20376135, -45147, -39828, 33990, -34509, 90294, 
  79332, -68304, 69342, 77, 68, -58, 59 ]);
subst := [&+[M[i,j]*R.j: j in [1..4]]:i in [1..4]];
N := -(1/(2^13*3^19))*Matrix(Rationals(),2,2,[173,-1/3,1773,-3]);
GG := [&+[N[i,j]*[G1,G2][j]^M: j in [1..2]]: i in [1..2]];
assert GG eq [ 
  -x^3 + 4*x^2*y + 3*x^2*z - x^2*t + 6*x*y^2 + 2*x*y*z - 2*x*y*t - 6*x*z^2 
  + 4*x*z*t - 11*x*t^2 + y^3 + 7*y^2*t - 2*y*z^2 + 4*y*z*t - 4*y*t^2 + 6*z^3 
  - 7*z^2*t + 4*z*t^2 + t^3, 2*x^3 - x^2*y + 5*x^2*t - 10*x*y^2 - 2*x*y*z 
  + 16*x*y*t - 3*x*z^2 + 4*x*z*t + 8*x*t^2 - 5*y^3 - y^2*z - 3*y^2*t 
  - y*z^2 - 2*y*z*t + 12*y*t^2 + 3*z^3 - 4*z^2*t + 2*z*t^2 - 3*t^3 ];
inv := Determinant(N)^6*Determinant(M)^9*inv;
assert inv eq -3^4*67^5;

P := [1,-2,-1,0];
c4 := -a/27;
c6 := -b/54;
assert Evaluate(GG[1],P) eq 0 and Evaluate(GG[2],P) eq 0;
P1 := [Evaluate(f,P): f in subst];
assert Evaluate(G1,P1) eq 0 and Evaluate(G2,P1) eq 0;
gamma1,gamma2 := Explode(TangentLineConstruction([G1,G2],P1));
la := gamma2;
mu := 3*gamma1;
DD,cc4,cc6 := HessePolynomials(3,1,[c4,c6]:Variables := [la,mu]);
c4star := -4*DD/(c4^3-c6^2);
c6star := -8*cc6/(c4^3-c6^2)^2;
E1 := MinimalModel(EllipticCurve([-27*c4star,-54*c6star]));
print CR(E1),aInvariants(E1);

print "Example 7.5";

QT<T> := FunctionField(Rationals());
a := (1/2)*(39*T^4 - 60*T^3 - 162*T^2 + 60*T + 39);
b := 47*T^6 + 120*T^5 + 21*T^4 + 21*T^2 - 120*T + 47;
RR<x,y,z,t> := PolynomialRing(QT,4);
F1,F2 := XE9Equations(a,b:R := RR);
P := [(15/2)*(3*T^4 + 8*T^3 - 2*T^2 - 8*T + 3),T^2 + 1,1,0];
assert Evaluate(F1,P) eq 0 and Evaluate(F2,P) eq 0;
gamma1,gamma2 := Explode(TangentLineConstruction([F1,F2],P));
la := 47*T^6 - 78*T^5 - 153*T^4 + 244*T^3 + 153*T^2 - 78*T - 47;
mu := 18*(T^2 + 1)*(T^2 + 6*T - 1);
assert la/mu eq gamma2/(3*gamma1);
c4 := -a/27;
c6 := -b/54;
DD,cc4,cc6 := HessePolynomials(3,1,[c4,c6]:Variables := [la,mu]);
E0 := EllipticCurve([Evaluate(a,0),Evaluate(b,0)]);
E1 := EllipticCurve([Evaluate(-27*cc4,0),Evaluate(-54*cc6,0)]);
print [Conductor(E0),Conductor(E1)];

print "Example 7.6";

QT<T> := FunctionField(Rationals());
a := 3*(3*T + 1)*(6*T^3 - 3*T - 1)*(9*T^3 - 9*T - 4)^2;
b := 2*(3*T^3 + 27*T^2 + 21*T + 4)*(6*T^3 - 3*T - 1)^2*(9*T^3 - 9*T - 4)^2; 
RR<x,y,z,t> := PolynomialRing(QT,4);
_,_,G1,G2 := XE9Equations(a,b:R := RR);
P := [-(6*T^3-3*T-1)*(9*T^3 - 9*T - 4),T,1,0];
assert Evaluate(G1,P) eq 0 and Evaluate(G2,P) eq 0;
gamma1,gamma2 := Explode(TangentLineConstruction([G1,G2],P));
la := (3*T+1)*(9*T^3-9*T-4)*(6*T^3-3*T-1)
                   *(180*T^4 + 321*T^3 + 216*T^2 + 66*T + 8);
mu := 3*(369*T^6 + 1107*T^5 + 1431*T^4 + 1017*T^3 + 414*T^2 + 90*T + 8);
assert la/mu eq gamma2/(3*gamma1);
c4 := -a/27;
c6 := -b/54;
DD,cc4,cc6 := HessePolynomials(3,1,[c4,c6]:Variables := [la,mu]);
c4star := -4*DD/(c4^3-c6^2);
c6star := -8*cc6/(c4^3-c6^2)^2;
E0 := EllipticCurve([a,b]);
E1 := EllipticCurve([-27*c4star,-54*c6star]);

Tlist := [-1,-1/3,-1/4,-1/2,-2/3];
dlist := [6,6,3,-6,6];
for i in [1..5] do
  EE0 := EllipticCurve([Evaluate(a,Tlist[i]): a in aInvariants(E0)]);
  EE0 := MinimalModel(QuadraticTwist(EE0,dlist[i]));
  EE1 := EllipticCurve([Evaluate(a,Tlist[i]): a in aInvariants(E1)]);
  EE1 := MinimalModel(QuadraticTwist(EE1,dlist[i]));
  print Tlist[i],dlist[i],CR(EE0),CR(EE1);
end for;

print "Example 7.7";

E := EllipticCurve("1782b1");
a := 765;
b := 15102;
assert IsIsomorphic(EllipticCurve([a,b]),E);
R<v,w,x,y,z> := PolynomialRing(Rationals(),5);
F := XE11Equations(a,b:R := R);
inv :=  -4*(4*a^3+27*b^2)^2;
assert inv eq -2^28*3^12*11^6;
mu := 2^19*3^4*11^5;
M := Matrix(Rationals(),5,5,[984, 12900, -9093, -34056, 13689,
  -2040, -24252, -3315, 0, -16857, 328, 164, -435, 0, -57, -352, 
  88, -264, 264, -1056, -8, -4, -13, 0, 25]);
subst := [&+[M[i,j]*R.j: j in [1..5]]:i in [1..5]];
FF0 := (1/mu)*(F^M);
assert FF0 eq -v^2*w + v^2*x - v^2*y + 2*v^2*z - v*w^2 + 4*v*w*z 
  - 4*v*x^2 - 8*v*x*y + 2*v*x*z + 6*v*y*z + 3*v*z^2 + 2*w^3 - 3*w^2*x
  - 2*w^2*y + 8*w^2*z + 6*w*x^2 + 2*w*x*y + 2*w*x*z + 6*w*y^2 - 6*w*y*z
  + 9*w*z^2 - x^3 - x^2*z - 3*x*y^2 - 6*x*y*z - 9*x*z^2 - 6*y^3 
  + 9*y^2*z + 3*y*z^2 - 7*z^3;
inv0 := (1/mu)^5*Determinant(M)^3*inv;
assert inv0 eq 2^2*3^4*11^2;

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

// load "quintics.m";

c4 := -a/27;
c6 := -b/54;
Delta := (c4^3 - c6^2)/1728;
FF,GG := XE11EquationsRaw(c4,Delta:R:=R);

v1 := c6*v + 2*c6*w - 6*c4^2*x + 3*c4^2*y - 9*c4*c6*z;
w1 := c6*v + 6*c4^2*x + 3*c4^2*y + 9*c4*c6*z;
v11 := 44*(2*c4*v - 6*c6*w + 33*c6*x + 135*c4^2*y + 810*c4*c6*z);
x11 := 60*(5*v + 729*c4*y + 2187*c6*z);
y11 := 11*(c4*v - 3*c6*w - 6*c6*x);
z11 := 60*(v + 27*c4*y + 81*c6*z);
subst1 := [-v1,w1,-864*x,-36*c4*x - 108*c6*z,72*y];
subst2 := [v11,-427680*y,x11,-y11,-z11];

M := Matrix(5,5,[MC(f,R.i): i in [1..5],f in subst1]);
assert Evaluate(FF,subst1) eq 8*c6^3*F;
inv1 := (8*c6^3)^5/Determinant(M)^3*inv;
cc4 := c4Invariant(FF);

H := Matrix(5,5,[Deriv(Deriv(F,i),j):i,j in [1..5]]);
HH := Matrix(5,5,[Deriv(Deriv(FF,i),j):i,j in [1..5]]);
HH0 := Matrix(5,5,[Deriv(Deriv(FF0,i),j):i,j in [1..5]]);

pts := [[-1,5,1,2,1],[0,0,0,1,0],[1,1,-1,0,-4]];
dd1 := d1(c4,Delta,[v,w,x,y,z]);
P4 := Proj(R);

for i in [1..#pts] do
  P0 := pts[i];
  assert forall{f : f in Minors(HH0,4) | Evaluate(f,P0) eq 0};
  P := [Evaluate(f,P0): f in subst];
  assert forall{f : f in Minors(H,4) | Evaluate(f,P) eq 0};
  P1 := [Evaluate(f,P): f in subst1];
  assert forall{f : f in Minors(HH,4) | Evaluate(f,P1) eq 0};
  FP := Evaluate(FF,P1);
  c4P := Evaluate(cc4,P1);
  j1 := c4P^3/(inv1^8*FP^11);
  dP := Evaluate(dd1,P1);
  if P4!P1 eq P4![1,0,0,0,0] then
    E1 := E;
    assert dP eq 0;
  else
    C4 := dP*c4P;
    DD := dP^3*inv1^8*FP^11;
    flag,C6 := IsSquare(C4^3 - 1728*DD);
    assert flag;
    E1 := MinimalModel(EllipticCurve([-27*C4,-54*C6]));
  end if;
  assert jInvariant(E1) eq j1;
  v := TracesOfFrobenius(E,100);
  v1 := TracesOfFrobenius(E1,100);
  mm := [GF(11)!(v[i]-v1[i]): i in [1..#v]];
  if #[m : m in mm | m ne 0 ] gt 5 then 
    E1 := MinimalModel(QuadraticTwist(E1,-1));
    v1 := TracesOfFrobenius(E1,100);
    mm := [GF(11)!(v[i]-v1[i]): i in [1..#v]];
  end if;
  print CR(E1),aInvariants(E1);
  print mm;
end for;

print "Example 7.8";

E := EllipticCurve("4466c1");
a := 85;
b := -83162;
assert IsIsomorphic(EllipticCurve([a,b]),E);
R<v,w,x,y,z> := PolynomialRing(Rationals(),5);
_,G := XE11Equations(a,b:R := R);
inv :=  8*(4*a^3+27*b^2);
assert inv eq 2^21*7*11^2*29^2;
mu := 2^23;
M := Matrix(Rationals(),5,5,[4096, -1408, 128, -1312, 45088, 0, 128, 128, 
  32, 110, 0, 0, -256, -96, -103, 0, 0, 0, -32, -11, 0, 0, 0, 0, -1]);
subst := [&+[M[i,j]*R.j: j in [1..5]]:i in [1..5]];
GG0 := (1/mu)*(G^M);
assert GG0 eq -2*v^2*z - 4*v*w*y + 12*v*x*y + 4*v*x*z + 5*v*y^2 + 6*v*y*z 
  - 43*v*z^2 - w^2*x + w^2*y - 4*w*x*y - 2*w*x*z - 3*w*y^2 + 196*w*y*z 
  + 83*w*z^2 - 11*x^3 - 12*x^2*y - 9*x^2*z - 11*x*y^2 + 366*x*y*z 
  + 125*x*z^2 + 322*y^3 + 447*y^2*z + 275*y*z^2 + 632*z^3;
inv0 := (1/mu)^5*Determinant(M)^3*inv;
assert inv0 eq -2^2*7*11^2*29^2;

c4 := -a/27;
c6 := -b/54;
Delta := (c4^3 - c6^2)/1728;
_,GG := XE11EquationsRaw(c4,Delta:R:=R);

v1 := c6*v + 2*c6*w - 6*c4^2*x + 3*c4^2*y - 9*c4*c6*z;
w1 := c6*v + 6*c4^2*x + 3*c4^2*y + 9*c4*c6*z;
v11 := 44*(2*c4*v - 6*c6*w + 33*c6*x + 135*c4^2*y + 810*c4*c6*z);
x11 := 60*(5*v + 729*c4*y + 2187*c6*z);
y11 := 11*(c4*v - 3*c6*w - 6*c6*x);
z11 := 60*(v + 27*c4*y + 81*c6*z);
subst1 := [-v1,w1,-864*x,-36*c4*x - 108*c6*z,72*y];
subst2 := [v11,-427680*y,x11,-y11,-z11];

M := Matrix(5,5,[MC(f,R.i): i in [1..5],f in subst2]);
assert Evaluate(GG,subst2) eq 2^5*3^6*5^3*11^3*c6^3*G;
inv1 := (2^5*3^6*5^3*11^3*c6^3)^5/Determinant(M)^3*inv;
cc4 := c4Invariant(GG);

H := Matrix(5,5,[Deriv(Deriv(G,i),j):i,j in [1..5]]);
HH := Matrix(5,5,[Deriv(Deriv(GG,i),j):i,j in [1..5]]);
HH0 := Matrix(5,5,[Deriv(Deriv(GG0,i),j):i,j in [1..5]]);

pts := [[-7,11,3,1,1],[7830,-3553,510,-281,71]];
dd2 := d2(c4,Delta,[v,w,x,y,z]);
P4 := Proj(R);

for i in [1..#pts] do
  P0 := pts[i];
  assert forall{f : f in Minors(HH0,4) | Evaluate(f,P0) eq 0};
  P := [Evaluate(f,P0): f in subst];
  assert forall{f : f in Minors(H,4) | Evaluate(f,P) eq 0};
  P1 := [Evaluate(f,P): f in subst2];
  assert forall{f : f in Minors(HH,4) | Evaluate(f,P1) eq 0};
  FP := Evaluate(GG,P1);
  c4P := Evaluate(cc4,P1);
  j1 := c4P^3/(inv1^8*FP^11);
  dP := Evaluate(dd2,P1);
  if dP eq 0 then 
    E1 := EllipticCurve("4466c2");
  else
    C4 := dP*c4P;
    DD := dP^3*inv1^8*FP^11;
    flag,C6 := IsSquare(C4^3 - 1728*DD);
    assert flag;
    E1 := MinimalModel(EllipticCurve([-27*C4,-54*C6]));
  end if;
  assert jInvariant(E1) eq j1;
  v := TracesOfFrobenius(E,100);
  v1 := TracesOfFrobenius(E1,100);
  mm := [GF(11)!(v[i]-v1[i]): i in [1..#v]];
  if #[m : m in mm | m ne 0 ] gt 5 then 
    E1 := MinimalModel(QuadraticTwist(E1,-1));
    v1 := TracesOfFrobenius(E1,100);
    mm := [GF(11)!(v[i]-v1[i]): i in [1..#v]];
  end if;
  print CR(E1),aInvariants(E1);
  print mm;
end for;

print "Example 7.9";

QT<T> := FunctionField(Rationals());

a := -3*(T - 3)*(T^4 - 5*T^2 - 24*T - 92)/(T^3 - T^2 + 4*T + 24);
b := -2*(T - 3)*(T^5 - T^4 - 11*T^3 - 43*T^2 
 		            - 62*T - 316)/(T^3 - T^2 + 4*T + 24);
A := -3*(T - 3)*(T^2 - 8*T - 17)*(T^3 - T^2 + 4*T + 24) 
  *(T^12 - 250*T^11 + 3473*T^10 - 23824*T^9 + 106654*T^8 - 354556*T^7 
    + 890186*T^6 - 1710568*T^5 + 2386357*T^4 - 2054170*T^3 + 1799781*T^2 
    + 956680*T + 3570796);
B := -2*(T - 3)*(T^3 - T^2 + 4*T + 24)^2*(T^20 + 476*T^19 - 27815*T^18 
  + 556718*T^17 - 6046664*T^16 + 42450848*T^15 - 213832636*T^14 
  + 823702888*T^13 - 2497998850*T^12 + 5954643736*T^11 
  - 10798748818*T^10  + 13644339892*T^9 - 7927895108*T^8 
  - 10398245632*T^7 + 25581636532*T^6 - 10366268760*T^5 - 60876061719*T^4 
  + 164062110060*T^3 - 98120800447*T^2 + 262948421518*T + 141270230564);

DeltaE := 2^12*3^6*(T - 5)*(T - 3)^2*(T + 1)^5
                           *(T^2 + 7)/(T^3 - T^2 + 4*T + 24)^3;
DeltaE1 := - 2^12*3^6*(T - 5)^4*(T - 3)^2*(T + 1)^3*(T^2 + 7)
               *(T^3 - T^2 + 4*T + 24)^3*(T^3 - T^2 + 15*T - 31)^11;

E := EllipticCurve([a,b]);
E1 := EllipticCurve([A,B]);
assert Discriminant(E) eq DeltaE;
assert Discriminant(E1) eq DeltaE1;

P := [T^6 + T^5 + 31*T^4 + 259*T^3 + 520*T^2 + 676*T + 1248,
  -(T - 3)*(T^5 + 4*T^4 + 43*T^3 + 100*T^2 - 44*T - 320),
  -(T^2 + 3*T + 14)*(T^3 - T^2 + 4*T + 24),
  0,
  (T + 4)*(T^3 - T^2 + 4*T + 24)];

RR<v,w,x,y,z> := PolynomialRing(QT,5);
F := XE11Equations(a,b:R := RR);
H := Matrix(5,5,[Deriv(Deriv(F,i),j):i,j in [1..5]]);
assert forall{f : f in Minors(H,4) | Evaluate(f,P) eq 0};

Tlist := [2,1,-3,11];
dlist := [-6,42,-2,-426];
for i in [1..4] do
  EE0 := EllipticCurve([Evaluate(a,Tlist[i]): a in aInvariants(E)]);
  EE0 := MinimalModel(QuadraticTwist(EE0,dlist[i]));
  EE1 := EllipticCurve([Evaluate(a,Tlist[i]): a in aInvariants(E1)]);
  EE1 := MinimalModel(QuadraticTwist(EE1,dlist[i]));
  print Tlist[i],dlist[i],CR(EE0),CR(EE1);
  v0 := TracesOfFrobenius(EE0,100);
  v1 := TracesOfFrobenius(EE1,100);
  print [GF(11)!(v0[i]-v1[i]): i in [1..#v0]];
end for;

