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

// VISIBLE 2-TORSION IN THE TATE-SHAFAREVICH GROUP OF AN ELLIPTIC CURVE

// Tom Fisher
// May 2026

// This is a Magma file accompanying the paper. It checks some of the
// calculations and examples, in the order they appear in the paper.

// Loading this file should produce no output -- meaning that
// all the "assert" statements are true.

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

// The invariants of a binary quartic

function Invs(BQ)
  a,b,c,d,e := Explode(Eltseq(BQ));
  I := 12*a*e - 3*b*d + c^2;
  J := 72*a*c*e - 27*a*d^2 - 27*b^2*e + 9*b*c*d - 2*c^3;
  return [I,J];
end function;

// Checking against the "GenusOneModel" functions

K<a,b,c,d,e,p,q,r,s> := FunctionField(Rationals(),9);
P<X> := PolynomialRing(K);
R2<x,z> := PolynomialRing(K,2);
f := a*x^4 + b*x^3*z + c*x^2*z^2 + d*x*z^3 + e*z^4;
BQ := GenusOneModel(f);
I,J := Explode(Invs(BQ));
assert cInvariants(BQ) eq [2^4*I, 2^5*J];
assert Discriminant(BQ) eq (16/27)*(4*I^3 - J^2);
E := EllipticCurve(X^3 - 27*I*X - 27*J);
assert IsIsomorphic(Jacobian(BQ),E);

// Checking against some older Magma functions

assert I eq QuarticIInvariant(Evaluate(f,[X, 1]));
assert J eq QuarticJInvariant(Evaluate(f,[X, 1]));

// Section 3 : Quadric intersections

Det := Determinant;
SymMat := func<Q|2*SymmetricMatrix(Q)>;

function AssociatedQuadraticForm(R3,BQ)
  a,b,c,d,e := Explode(Eltseq(BQ));
  x0,x1,x2 := Explode([R3.i: i in [1..3]]);
  return a*x0^2 + b*x0*x1 + (c/3)*(x0*x2 + 2*x1^2) + d*x1*x2 + e*x2^2;
end function;

R3<x0,x1,x2> := PolynomialRing(K,3);
Q := x0*x2 - x1^2;
Qf := AssociatedQuadraticForm(R3,BQ);
assert Evaluate(Qf,[x^2, x*z, z^2]) eq f;

// Checking the proof of Lemma 3.3

assert Det(X*SymMat(Q) - 9*SymMat(Qf)) eq 2*(X^3 - 27*I*X - 27*J);
f1 := Evaluate(f,[p*x + r*z, q*x + s*z]);
BQ1 := GenusOneModel(f1);
Qf1 := AssociatedQuadraticForm(R3,BQ1);
xx0 := p^2*x0 + 2*p*r*x1 + r^2*x2;
xx1 := p*q*x0 + (p*s + q*r)*x1 + r*s*x2;
xx2 := q^2*x0 + 2*q*s*x1 + s^2*x2;
assert Evaluate(Qf,[xx0, xx1, xx2]) eq Qf1;
assert Evaluate(Q,[xx0, xx1, xx2]) eq (p*s - q*r)^2*Q;

// Definition 3.4

function AssociatedQuadricIntersection(R6,f,g)
  u0,u1,u2,v0,v1,v2 := Explode([R6.i : i in [1..6]]);
  K := BaseRing(f);
  R3 := PolynomialRing(K,3);
  Qf := AssociatedQuadraticForm(R3,f);
  Qg := AssociatedQuadraticForm(R3,g);
  F := u0*u2 - u1^2 - v0*v2 + v1^2;
  G := Evaluate(Qf,[u0, u1, u2]) - Evaluate(Qg,[v0, v1, v2]);
  return [F, G];
end function;

// Checking det(x*F - 9*G) = -4*(X^3 - 27*I*X - 27*J)^2

K<a,b,c,d,e,aa,bb,cc,dd,ee> := FunctionField(Rationals(),10);
R2<x,z> := PolynomialRing(K,2);
P<X> := PolynomialRing(K);
f := GenusOneModel(a*x^4 + b*x^3*z + c*x^2*z^2 + d*x*z^3 + e*z^4);
g := GenusOneModel(aa*x^4 + bb*x^3*z + cc*x^2*z^2 + dd*x*z^3 + ee*z^4);
I,J := Explode(Invs(f));
II,JJ := Explode(Invs(g));
R6<u0,u1,u2,v0,v1,v2> := PolynomialRing(K,6);
F,G := Explode(AssociatedQuadricIntersection(R6,f,g));
detform := Determinant(X*SymMat(F) - 9*SymMat(G));
assert detform eq -4*(X^3 - 27*I*X - 27*J)*(X^3 - 27*II*X - 27*JJ);

// Checking Remark 3.5

// Over an algebraically closed field we may assume that f = g and
// that this binary quartic is in "Legendre form".

K<t> := FunctionField(Rationals());
R2<x,z> := PolynomialRing(K,2);
R6<x1,x2,x3,x4,x5,x6> := PolynomialRing(K,6);
f := GenusOneModel(x*z*(x - z)*(x - t*z));
eqns := AssociatedQuadricIntersection(R6,f,f);

// A change of coordinates that rewrites our quadric intersection
// in the required form. We may assume that 6*t*(t - 1) <> 0.

A := Matrix(K,6,6,[
   [ t, t, t^2, t^2, t^2 - t, t^2 - t ],
   [ t, t, t, t, 0, 0 ],
   [ 1, 1, t, t, -t + 1, -t + 1 ],
   [ t, -t, t^2, -t^2, t^2 - t, -t^2 + t ],
   [ t, -t, t, -t, 0, 0 ],
   [ 1, -1, t, -t, -t + 1, t - 1 ] ]);
B := Matrix(K,2,2,[ t^2 - 1, 3*t - 3, -2*t + 1, 3 ]);
assert Determinant(A) eq 2^5*t^4*(t - 1)^4;
assert Determinant(B) eq 3^2*t*(t - 1);
B *:= -1/(12*t^2*(t - 1)^2);
eqns := [q^A: q in eqns];
eqns := [&+[B[i,j]*eqns[j]: j in [1..2]]: i in [1..2]];
assert eqns eq [ x1*x2 - x3*x4, x3*x4 - x5*x6 ];

// Let X = { x1*x2 = x3*x4 = x5*x6 } subset P^5. We claim that
// the only planes contained in X are the obvious ones, i.e.,
// P_ijk = {x_i = x_j = x_k = 0} for i in {1,2}, j in {3,4}, k in {5,6}.

R6<x1,x2,x3,x4,x5,x6> := PolynomialRing(Rationals(),6);
X := Scheme(Proj(R6),[ x1*x2 - x3*x4, x3*x4 - x5*x6 ]);
P_235 := Scheme(X,[x2, x3, x5]);
P_236 := Scheme(X,[x2, x3, x6]);
P_245 := Scheme(X,[x2, x4, x5]);
P_246 := Scheme(X,[x2, x4, x6]);
line := Scheme(X,[x1 - x3, x3 - x5, x2 - x4, x4 - x6]);

// Let Lambda be a plane contained in X. It must meet {x4 = x6 = 0},
// and so (up to permuting and rescaling the coordinates) contains
// one of the following points P.

P := X![1, 0, 1, 0, 1, 0];
assert (X meet TangentSpace(P)) eq (line join P_246);
P := X![1, 0, 1, 0, 0, 0];
assert (X meet TangentSpace(P)) eq (P_245 join P_246);
P := X![1, 0, 0, 0, 0, 0];
assert (X meet TangentSpace(P)) eq &join[P_235,P_236,P_245,P_246];

// Since Lambda must be contained in the tangent space at P,
// this proves the claim.

// The planes arising from proper equivalences

K<a,b,c,d,e,p,q,r,s> := FunctionField(Rationals(),9);
R2<x,z> := PolynomialRing(K,2);
g := a*x^4 + b*x^3*z + c*x^2*z^2 + d*x*z^3 + e*z^4;
f := Evaluate(g,[p*x + r*z,q*x + s*z])/(p*s - q*r)^2;
f := GenusOneModel(f);
g := GenusOneModel(g);
R6<u0,u1,u2,v0,v1,v2> := PolynomialRing(K,6);
X_fg := Scheme(Proj(R6),AssociatedQuadricIntersection(R6,f,g));
Lambda := Scheme(Ambient(X_fg),[
  (p*s - q*r)*v0 - (p^2*u0 + 2*p*r*u1 + r^2*u2),
  (p*s - q*r)*v1 - (p*q*u0 + (p*s + q*r)*u1 + r*s*u2),
  (p*s - q*r)*v2 - (q^2*u0 + 2*q*s*u1 + s^2*u2)]);
assert Lambda subset X_fg;

// Checking the proof of Lemma 3.6

R<p,q,r> := PolynomialRing(Rationals(),3);
P<X> := PolynomialRing(R);
s := -p; // trace = 0
det := p*s - q*r;
M := Matrix(3,3,[
  [ p^2, 2*p*r, r^2 ],
  [ p*q, -p^2 + q*r, -p*r ],
  [ q^2, -2*p*q, p^2 ] ]);
assert P!CharacteristicPolynomial(M) eq (X - det)*(X + det)^2;
N := Matrix(3,4,[
  [ -r, 0, r, 2*p ],
  [ p, r, 0, q ],
  [ q, -2*p, q, 0 ]]);
assert GCD(Minors(N,3)) eq -det;
assert M*N eq det*N*DiagonalMatrix(R,[1,-1,-1,-1]);

// The first column of N is the only eigenvector of M with
// eigenvalue p*s - q*r. This shows that the planes corresponding
// to [1,0;0,1] and [p,q;r,s] meet only at (-r:p:q:-r:p:q).

// Next we check that this is a singular point.
// (This could also be deduced from Remark 3.5.)

K<a,b,c,d,e,p,q,r> := PolynomialRing(Rationals(),8);
s := -p; // trace = 0
R2<x,z> := PolynomialRing(K,2);
f := a*x^4 + b*x^3*z + c*x^2*z^2 + d*x*z^3 + e*z^4;
f1 := Evaluate(f,[p*x + r*z,q*x + s*z]);
relns := Coefficients(f1 - (p*s - q*r)^2*f);
R3<x0,x1,x2> := PolynomialRing(K,3);
Q := x0*x2 - x1^2;
Qf := AssociatedQuadraticForm(R3,GenusOneModel(f));
Jmat := Matrix(2,3,[Derivative(q,j): j in [1..3],q in [Q,Qf]]);
Jmat := Evaluate(Jmat,[-r,p,q]);
M := Matrix(3,5,[
  [ 0, 6*r^2, -8*p*r, 6*p^2, 0 ],
  [ 12*r^2, -3*p*r, 2*q*r, -3*p*q, 0 ],
  [ 24*p*r, -6*p^2 + 6*q*r, -4*p*q, 0, 0 ] ]);
assert [6*(p*s - q*r)^2*m : m in Minors(Jmat,2)] eq
  [&+[M[i,j]*relns[j] : j in [1..5]]: i in [1..3]];

// Checking the proof of Lemma 3.7

K<la,mu> := FunctionField(Rationals(),2);
R6<z12,z13,z14,z23,z24,z34> := PolynomialRing(K,6);
Z := Matrix(4,4,[0,z12,z13,z14,0,0,z23,z24,0,0,0,z34,0,0,0,0]);
Z -:= Transpose(Z);
F := z12*z34 - z13*z24 + z14*z23;
assert F eq Pfaffian(Z);
G := la*z12*z34 + mu*z13*z24;
R<x1,x2,x3,x4> := PolynomialRing(K,4);
RR<y1,y2,y3,y4> := PolynomialRing(R,4);
x := [x1,x2,x3,x4];
y := [y1,y2,y3,y4];
subst := [x[i]*y[j] - x[j]*y[i] where i,j is Explode(a): a in aa]
  where aa is [[1, 2], [1, 3], [1, 4], [2, 3], [2, 4], [3, 4]];
mat := SymMat(Evaluate(G,subst));
assert mat eq Matrix(R,4,4,[
  [ 0, mu*x3*x4, la*x2*x4, -(la + mu)*x2*x3 ],
  [ mu*x3*x4, 0, -(la + mu)*x1*x4, la*x1*x3 ],
  [ la*x2*x4, -(la + mu)*x1*x4, 0, mu*x1*x2 ],
  [ -(la + mu)*x2*x3, la*x1*x3, mu*x1*x2, 0 ]]);
assert Vector([x1, x2, x3, x4])*mat eq Vector([0, 0, 0, 0]);
assert GCD(Minors(mat,3)) eq x1*x2*x3*x4;

// Checking the proof of Proposition 5.1

R6<x12,x13,x14,x23,x24,x34> := PolynomialRing(Rationals(),6);
eqns := [x12*x34 - x13*x24,x13*x24 - x14*x23];
X := Scheme(Proj(R6),eqns);
Lambda1 := Scheme(X,[x12,x13,x14]);
Lambda2 := Scheme(X,[x12,x23,x24]);
Lambda3 := Scheme(X,[x13,x23,x34]);
Lambda4 := Scheme(X,[x14,x24,x34]);
mb := MinimalBasis(&join[Lambda1,Lambda2,Lambda3,Lambda4]);
assert {* Degree(f): f in mb *} eq {* 2^^3, 3^^4 *};
aa := [ x12*x34, x13*x24, x14*x23 ];
M := [Matrix(R6,4,4,[
  [ x12*x13*x14, alpha*x12, alpha*x13, alpha*x14 ],
  [ alpha*x12, x12*x23*x24, alpha*x23, alpha*x24 ],
  [ alpha*x13, alpha*x23, x13*x23*x34, alpha*x34 ],
  [ alpha*x14, alpha*x24, alpha*x34, x14*x24*x34 ]])
  : alpha in aa];
bb := [ M[1,i,i] : i in [1..4]];
assert Set(mb) eq Set(aa) join Set(bb);
cubics1 := [ a*R6.j : a in aa,j in [1..6]] cat bb; // vanishing on Lambdas
cubics2 := [ q*R6.j : q in eqns, j in [1..6]];     // vanishing on X
cubics3 := [ M[1,i,j] : i,j in [1..4] | i le j];   // complement
mons := MonomialsOfDegree(R6,3);
V := VectorSpace(Rationals(),#mons);
U1 := sub<V|[V![MonomialCoefficient(f,m): m in mons]: f in cubics1]>;
U2 := sub<V|[V![MonomialCoefficient(f,m): m in mons]: f in cubics2]>;
U3 := sub<V|[V![MonomialCoefficient(f,m): m in mons]: f in cubics3]>;
assert Dimension(U1) eq 22;
assert Dimension(U2) eq 12;
assert Dimension(U3) eq 10;
assert U1 eq U2 + U3;
assert Dimension(U2 meet U3) eq 0;
P3<x1,x2,x3,x4> := ProjectiveSpace(Rationals(),3);
R10<x11,x12,x13,x14,x22,x23,x24,x33,x34,x44>
  := PolynomialRing(Rationals(),10);
MM := Matrix(R10,4,4,[
  [x11, x12, x13, x14],
  [x12, x22, x23, x24],
  [x13, x23, x33, x34],
  [x14, x24, x34, x44]]);
Y := Scheme(Proj(R10),Minors(MM,2));
phi := map<X->Y|[M[1,i,j] : i,j in [1..4] | i le j]>;
psi := map<P3->X|[x1*x2, x1*x3, x1*x4, x2*x3, x2*x4, x3*x4]>;
iota := iso<P3 -> Y|[P3.i*P3.j: i,j in [1..4] | i le j],[x11, x12, x13, x14]>;
assert psi*phi eq iota; // beware Magma's order convention

// Checking phi is regular at all smooth points.

R6<x12,x13,x14,x23,x24,x34> := R6;
N := [x13*x14*M[1],x12*x14*M[2],x12*x13*M[3],x23*x24*M[1]];
N[1,3,3] +:= x13^2*x34*(eqns[1] + eqns[2]);
N[1,4,4] +:= x14^2*x34*eqns[1];
N[2,2,2] +:= x12^2*x24*eqns[2];
N[2,4,4] -:= x14^2*x24*eqns[1];
N[3,2,2] -:= x12^2*x23*eqns[2];
N[3,3,3] -:= x13^2*x23*(eqns[1] + eqns[2]);
N[4,3,3] +:= x23^2*x34*eqns[1];
N[4,4,4] +:= x24^2*x34*(eqns[1] + eqns[2]);
N := [M[1]] cat [ChangeRing(1/[x12,x13,x14,x12][i]*N[i],R6) : i in [1..4]];
phi := map<X->Y|[[N[k,i,j] : i,j in [1..4] | i le j] : k in [1..5]]>;
assert ReducedSubscheme(BaseScheme(phi)) eq SingularSubscheme(X);

// Some details in the proof of Proposition 6.2.

K<a,b,c,d,e,la,mu,la1,mu1> := FunctionField(Rationals(),9);
P<X> := PolynomialRing(K);
R2<x,z> := PolynomialRing(K,2);
R3 := PolynomialRing(K,3);
R4<x0,x1,x2,y> := PolynomialRing(K,4);
det := la*mu1 - mu*la1;
BQ := GenusOneModel(a*x^4 + b*x^3*z + c*x^2*z^2 + d*x*z^3 + e*z^4);
Qf := AssociatedQuadraticForm(R3,BQ);
I,J := Explode(Invs(BQ));
F := mu*(x0*x2 - x1^2) - la*Evaluate(Qf,[x0,x1,x2]);
G := mu1*(x0*x2 - x1^2) - la1*Evaluate(Qf,[x0,x1,x2]) - y^2;
f := X^3 - 27*I*X - 27*J;
f1 := -(9*mu*X + 9*mu1)^3 + 27*I*(9*mu*X + 9*mu1)*(la*X + la1)^2
  + 27*J*(la*X + la1)^3;
assert Determinant(X*SymMat(F) + SymMat(G)) eq (2/27)^2*f1;
X1 := -(la1*X - 9*mu1)/(la*X - 9*mu);
assert Evaluate(f1,X1) eq -3^6*det^3*f/(la*X - 9*mu)^3;

// Example 7.1

R<x,y> := PolynomialRing(Rationals(),2);
Weqn := func<E|Evaluate(Equation(E),[x, y, 1])>;

E := EllipticCurve("571a1");
assert Weqn(E) eq (y^2 + y) - (x^3 - x^2 - 929*x - 10595);
assert Rank(E) eq 0;
assert #TorsionSubgroup(E) eq 1;
analytic_order_sha := ConjecturalRegulator(E);
assert Abs(analytic_order_sha - 4) lt 10^(-20);

F := EllipticCurve("571b1");
assert Weqn(F) eq (y^2 + y) - (x^3 + x^2 - 4*x + 2);
assert Rank(F) eq 2;
assert #TorsionSubgroup(F) eq 1;

R2<x,z> := PolynomialRing(Rationals(),2);
P<X> := PolynomialRing(Rationals());
f := -11*x^4 + 68*x^3*z - 52*x^2*z^2 - 164*x*z^3 - 64*z^4;
g := -4*x^4 - 60*x^3*z - 232*x^2*z^2 - 52*x*z^3 - 3*z^4;
h := -31*x^4 - 78*x^3*z + 32*x^2*z^2 + 102*x*z^3 - 53*z^4;
BQlist := [GenusOneModel(m): m in [f,g,h]];
I := 44608;
J := 18842960;
assert IsIsomorphic(E,EllipticCurve(X^3 - 27*I*X - 27*J));
assert forall{BQ : BQ in BQlist | Invs(BQ) eq [I,J]};

// Checking agreement with 2-descent

td := TwoDescent(E);
BQlist1 := [GenusOneModel(C): C in td];
mat := Matrix(3,3,[IsEquivalent(BQ,BQ1) select 1 else 0:
   BQ1 in BQlist1,BQ in BQlist]);
assert Vector([1, 1, 1])*mat eq Vector([1, 1, 1]);
assert Vector([1, 1, 1])*Transpose(mat) eq Vector([1, 1, 1]);

// A rational point on X_{f,g}

R6<u0,u1,u2,v0,v1,v2> := PolynomialRing(Rationals(),6);
P5 := Proj(R6);
eqns := AssociatedQuadricIntersection(R6,BQlist[1],BQlist[2]);
X_fg := Scheme(P5,eqns);
pt := X_fg![ 3, 1, 0, -1, 0, 1 ];
subst := [pt[i] : i in [1..3]] cat [0, 0, 0];
la := Evaluate(eqns[1],subst);
mu := Evaluate(eqns[2],subst);
la1 := 0;
mu1 := -1;
assert [la, mu, la1, mu1] eq [-1, 211/3, 0, -1];
assert la*mu1 - mu*la1 eq 1;
f1 := -(9*mu*X + 9*mu1)^3 + 27*I*(9*mu*X + 9*mu1)*(la*X + la1)^2
  + 27*J*(la*X + la1)^3;
c := Coefficient(f1,3);
F1 := EllipticCurve(c^2*Evaluate(f1,X/c));
assert IsIsomorphic(F,F1);

// Example 7.2

E := EllipticCurve("192a2");
R<x,y> := R;
assert Weqn(E) eq y^2 - (x - 1)*(x - 3)*(x + 3);
assert Rank(E) eq 1;
assert #TorsionSubgroup(E) eq 4;
R2<x,z> := R2;
f := -x^4 - 4*x^2*z^2 - z^4;
g := x^4 + 2*x^3*z + 2*x^2*z^2 + 4*x*z^3 + 4*z^4;
BQlist := [GenusOneModel(m): m in [f,g]];
I := 28;
J := -160;
assert 4*I^3 - J^2 eq 2^8*3^5;
assert IsIsomorphic(E,EllipticCurve(X^3 - 27*I*X - 27*J));
assert forall{BQ : BQ in BQlist | Invs(BQ) eq [I,J]};
eqns := AssociatedQuadricIntersection(R6,BQlist[1],BQlist[2]);
eqns := [X*eqns[1] - eqns[2] : X in [2/3, -10/3]];
assert eqns eq [
  (u0 + u2)^2 + 2*u1^2 + (v0 + v1)^2 + (v1 + 2*v2)^2,
  (u0 - u2)^2 + 6*u1^2 + (v0 + v1 + 2*v2)^2 - 3*v1^2 ];
XX := Scheme(P5,eqns);
SS := SingularSubscheme(XX);
assert Dimension(SS) eq 0;
assert Degree(SS) eq 6;
K<rt2> := QuadraticField(2);
pt1 := SS(K)![rt2, 0, -rt2, 2, -2, 1];
K<rtm2> := QuadraticField(-2);
pt2 := SS(K)![rtm2, 0, rtm2, -2, 0, 1];
K<i> := QuadraticField(-1);
pt3 := SS(K)![0, i, 0, 2, -1, 1];
K<alpha> := NumberField(X^4 + X^3 - X^2 + X + 4);
pt4 := XX(K)![alpha^2 - 1, 1, 0, alpha^2 - 1, alpha, 1];
assert not IsSingular(pt4);

// Example 7.3

K<p> := PolynomialRing(Rationals());
R2<x,z> := PolynomialRing(K,2);
f := 20*x^2*z^2 + p*(x^4 + 30*x^2*z^2 + 5*z^4);
g := (x^2 - 5*z^2)^2 + p*(x^4 + 4*x^3*z - 30*x^2*z^2 + 20*x*z^3 + 25*z^4);
BQlist := [GenusOneModel(m): m in [f,g]];
I := 2^4*5*(12*p^2 + 15*p + 5);
J := -2^6*5^2*(3*p + 2)*(9*p^2 + 15*p + 5);
assert 4*I^3 - J^2 eq 2^12*3^3*5^3*p^2*(11*p^2 + 15*p + 5)^2;
assert forall{BQ : BQ in BQlist | Invs(BQ) eq [I,J]};

R6<u0,u1,u2,v0,v1,v2> := PolynomialRing(K,6);
eqns := AssociatedQuadricIntersection(R6,BQlist[1],BQlist[2]);
eqns[2] := eqns[2] - 20/3*eqns[1];
H0 := 20*u1^2 - (v0 - 5*v2)^2;
H1 := u0^2 + 10*u0*u2 + 20*u1^2 + 5*u2^2;
H2 := v0^2 + 4*v0*v1 - 10*v0*v2 - 20*v1^2 + 20*v1*v2 + 25*v2^2;
assert eqns[2] eq H0 + p*(H1 - H2);
R6<x1,x2,x3,x4,x5,x6> := R6;
subst := [x1, p*x5, x2, 5*x4 + p*x6, x3, x4];
eqns := [Evaluate(q,subst): q in eqns];
eqns[2] := ExactQuotient(eqns[2],p);
F := x1*x2 + x3^2 - 5*x4^2;
G := x1^2 + 10*x1*x2 + 5*x2^2 + 20*x3^2 - 40*x3*x4;
F1 := -x4*x6 - p*x5^2;
G1 := -4*x3*x6 + 20*(p + 1)*x5^2 - (p + 1)*x6^2;
assert eqns eq [F + p*F1,G + p*G1];
K<z> := CyclotomicField(5);
rt5 := 1 + 2*(z + z^4);
assert rt5^2 eq 5;
RR<x1,x2,x3,x4,x5,x6> := PolynomialRing(K,6);
assert (-10 + 2*rt5)*(RR!F) + (RR!G) eq
  (x1 + rt5*x2)^2 - (z - z^4)^2*(2*x3 - (5 - rt5)*x4)^2;

// Some details in the proof of Proposition 8.1.

K<a11,a12,a13,a21,a22,a23,a31,a32,a33> := FunctionField(Rationals(),9);
P1P1<x1,x2,y1,y2> := ProductProjectiveSpace(K,[1,1]);
A := Matrix(K,3,3,[K.i : i in [1..9]]);
AA := ChangeRing(A,CoordinateRing(P1P1));
F := (Matrix(1,3,[x1^2,x1*x2,x2^2])*AA*Matrix(3,1,[y1^2,y1*y2,y2^2]))[1,1];

R<x,z> := PolynomialRing(K,2);
P<X> := PolynomialRing(R);
f := GenusOneModel(Discriminant(Evaluate(F,[x,z,X,1])));
g := GenusOneModel(Discriminant(Evaluate(F,[X,1,x,z])));
R6<u0,u1,u2,v0,v1,v2> := PolynomialRing(K,6);
eqns := AssociatedQuadricIntersection(R6,f,g);
l := Eltseq(Matrix(1,3,[u0,u1,u2])*ChangeRing(A,R6));
m := Eltseq(ChangeRing(A,R6)*Matrix(3,1,[v0,v1,v2]));
kappa := 8*a11*a33 - 4*a12*a32 + 8*a13*a31 - 4*a21*a23 + a22^2;
eqns[2] -:= (kappa/3)*eqns[1];
assert eqns eq [
  u0*u2 - u1^2 - (v0*v2 - v1^2),
  l[2]^2 - 4*l[1]*l[3] - (m[2]^2 - 4*m[1]*m[3])];

P3<x,y,z,t> := ProjectiveSpace(K,3);
uu0 := -a21*x^2 - a22*x*y - a23*y^2 - 2*a31*x*z - a32*(x*t + y*z) - 2*a33*y*t;
uu1 := a11*x^2 + a12*x*y + a13*y^2 - a31*z^2 - a32*z*t - a33*t^2;
uu2 := a21*z^2 + a22*z*t + a23*t^2 + 2*a11*x*z + a12*(x*t + y*z) + 2*a13*y*t;
vv0 := -a12*x^2 - a22*x*z - a32*z^2 - 2*a13*x*y - a23*(x*t + y*z) - 2*a33*z*t;
vv1 := a11*x^2 + a21*x*z + a31*z^2 - a13*y^2 - a23*y*t - a33*t^2;
vv2 := a12*y^2 + a22*y*t + a32*t^2 + 2*a11*x*y + a21*(x*t + y*z) + 2*a31*z*t;
param := [uu0,uu1,uu2,vv0,vv1,vv2];
assert forall{q : q in eqns | Evaluate(q,param) eq 0};
X := Scheme(Proj(R6),eqns);
psi := map<P3->X|param>;
assert Dimension(BaseScheme(psi)) eq 0;
assert Degree(BaseScheme(psi)) eq 4;

// END OF FILE