// Magma functions for the paper // "Diagonal genus 5 curves, elliptic curves over \Q(t), and rational diophantine quintuples" // by Michael Stoll // // This file contains the functions that implement the algorithm described // at the end of Section 2, for curves arising from (rational) diophantine quadruples. // To speed up (or make feasible) the class group computations // for the octic fields coming up in the Selmer group computations, // we allow Magma to assume GRH. // (See diophtuples-verify.magma for functionality to verify the // results without this assumption.) SetClassGroupBounds("GRH"); printf "\n======== USING GRH! =========\n\n"; // We need the polynomial ring over Q. P := PolynomialRing(Rationals()); // Determine models of the five elliptic curves arising as quotients // of the diagonal genus 5 curve associated to the diophantine // quadruple (a_1, a_2, a_3, a_4). function ellipticcurves(aseq) // aseq = [a1, a2, a3, a4] a,b,c,d := Explode(aseq); // set up P^3 as ambient space for intersections of two quadrics Pr3 := ProjectiveSpace(Universe(aseq), 3); // u1, ..., u4 are the projective coordinates u1,u2,u3,u4 := Explode([Pr3.i : i in [1..4]]); // eliminating one of the five coordinates from the three quadrics // defining the (projective) genus 5 curve gives two quadrics in four vars Cseq := [Curve(Pr3, [b*u2^2-c*u1^2-(b-c)*u4^2, b*u3^2-d*u1^2-(b-d)*u4^2]), Curve(Pr3, [a*u2^2-c*u1^2-(a-c)*u4^2, a*u3^2-d*u1^2-(a-d)*u4^2]), Curve(Pr3, [a*u2^2-b*u1^2-(a-b)*u4^2, a*u3^2-d*u1^2-(a-d)*u4^2]), Curve(Pr3, [a*u2^2-b*u1^2-(a-b)*u4^2, a*u3^2-c*u1^2-(a-c)*u4^2]), Curve(Pr3, [(c-a)*u2^2-(c-b)*u1^2-(b-a)*u3^2, (d-a)*u2^2-(d-b)*u1^2-(b-a)*u4^2])]; // collect results in a Magma list resseq := [**]; aaseq := Append(aseq, 0); // "homogenized" quadruple for i := 1 to #Cseq do // get isomorphisms with (standard) elliptic curves mapping (1:1:1:1) to the origin C := Cseq[i]; E1, toE1 := EllipticCurve(C, C![1,1,1,1]); // move to elliptic curve with explicit and simple equation as := [aaseq[j] : j in [1..5] | j ne i]; part := [as[1]*as[2] + as[3]*as[4], as[1]*as[3] + as[2]*as[4], as[1]*as[4] + as[2]*as[3]]; E := EllipticCurve(HyperellipticCurve(&*[Polynomial([p,1]) : p in part])); flag, toE := IsIsomorphic(E1, E); assert flag; // the entries of resseq are triples , where // C is the quadric intersection of genus 1, // E is the elliptic curve and // CtoE is the isomorphism C --> E Append(~resseq, ); end for; return resseq; end function; // Determine the squarefree integer in the square class of a rational number function sqfp(n) s := SquarefreeFactorization(Numerator(n)*Denominator(n)); return s; end function; // Define P^1 over Q (neede for some computations later) Pr1Q := ProjectiveSpace(Rationals(), 1); // Construct the projective diagonal genus 2 curve associated // to the diophantine quadruple (a, b, c, d) // (the additional "homogenizing" variable is last). function genus5curve(aseq) a,b,c,d := Explode(aseq); Pr4 := ProjectiveSpace(Universe(aseq), 4); return Curve(Pr4, [b*u1^2-a*u2^2+(a-b)*u0^2, c*u1^2-a*u3^2+(a-c)*u0^2, d*u1^2-a*u4^2+(a-d)*u0^2]); end function; // Construct the elliptic curves F_{T,\xi} over degree 4 fields // defined in Section 2, where \xi is the class of the point pt. // This is invariant under scaling the entries of the quadruple, // so we can assume that they are coprime integers. function ellcdeg4(aseqn, pt) // aseqn = [a_1,a_2,a_3,a_4,0] coprime integers // pt = point on C determining the relevant twists u := [R.i : i in [1..5]] where R := CoordinateRing(Ambient(Curve(pt))); function param(i,j,k) // give parameterization of the conic // (a_j-a_k) u_i^2 + (a_k-a_i) u_j^2 + (a_i-a_j) u_k^2 = 0 // in terms of dehomogenized binary quadratic forms a := aseqn[j]-aseqn[k]; b := aseqn[k]-aseqn[i]; c := aseqn[i]-aseqn[j]; // this can be done explicitly: par := [b*(x+1)^2 + c*x^2, c + a*(x+1)^2, a*x^2 + b]; assert a*par[1]^2 + b*par[2]^2 + c*par[3]^2 eq 0; // the second return value is x in terms of the projective coordinates return par, b/a * (u[j]+u[k])/(u[i]+u[k]); end function; function rad(i,j,k) // give coefficient of -u_k^2 in *u_i^2 = u_j^2 - *u_k^2 return (aseqn[j]-aseqn[i])/(aseqn[k]-aseqn[i]); end function; // we collect the results in a Magma list Eseq Eseq := [**]; // info is a parallel sequence that stores the combinatorial information // (the index pairs giving D_ij, D_ik in Case 2, the pairs giving D_ij, D_kl // together with {i,k} in Case 1) info := []; // loop over all unordered pairs of index pairs for pair in Subsets(Subsets({1..5}, 2), 2) do pair1, pair2 := Explode(Setseq(pair)); if IsEmpty(pair1 meet pair2) then // case 1: pair1 = {i,j}, pair2 = {k,l}, m is remaining index m := Rep({1..5} diff (pair1 join pair2)); // run through the four possible choices of i and k for i in pair1 do j := Rep(pair1 diff {i}); for k in pair2 do l := Rep(pair2 diff {k}); // parameterize the first conic par, xcoord := param(j,l,m); par := [p/g : p in par] where g := GCD(&cat[ChangeUniverse(Coefficients(p), Integers()) : p in par]); // determine the biquadratic field K rad1 := rad(k,m,l); rad2 := rad(i,m,j); r1 := sqfp(rad1); r2 := sqfp(rad2); if r1 eq 1 or r2 eq 1 or r1 eq r2 then // K is not of degree 4; skip // (instead of throwing these away, we should work with the two // curves over quadratic fields or four curves over Q, but this // would complicate the data structures and the program flow) printf "ellcdeg4: %o, %o --> smaller degree fields (sqrt(%o), sqrt(%o))\n", pair1, pair2, r1, r2; continue; end if; // set up K and the polynomial ring over K K := AbsoluteField(ext); X := PolynomialRing(K).1; // find the square roots of b and b' as elements of K sqrt1 := Roots(X^2 - rad1)[1,1]; sqrt2 := Roots(X^2 - rad2)[1,1]; // plug the parameterization into the second conic parK := [Evaluate(p, X) : p in par]; pol := Denominator(s1)*Denominator(s2)*(parK[3] - sqrt1*parK[2])*(parK[3] - sqrt2*parK[1]) where _, s1 := IsSquare(rad1/r1) where _, s2 := IsSquare(rad2/r2); // get the correct twist ptx := Evaluate(xcoord, Eltseq(pt)); polpt := Evaluate(pol, ptx); // construct F_{T,\xi} as a genus 1 "hyperelliptic" curve HE := HyperellipticCurve(polpt*pol); // map points from the genus 5 curve (actually, this gives the sequence of points on HE // with the same x-coordinate as the image) toHE := func; // the x-coordinate map on HE HEcov := map Pr1Q | [HE.1, HE.3]>; Append(~Eseq, ); Append(~info, [pair1, pair2, {i, k}]); end for; end for; else // case 2: pair1 = {i,j}, pair2 = {i,k}, l,m the remaining indices i := Rep(pair1 meet pair2); j := Rep(pair1 diff {i}); k := Rep(pair2 diff {i}); l, m := Explode(Setseq({1..5} diff (pair1 join pair2))); // parameterize the first conic par, xcoord := param(i,l,m); par := [p/g : p in par] where g := GCD(&cat[ChangeUniverse(Coefficients(p), Integers()) : p in par]); // determine K etc. (see case 1) rad1 := rad(j,l,m); rad2 := rad(k,l,m); r1 := sqfp(rad1); r2 := sqfp(rad2); if r1 eq 1 or r2 eq 1 or r1 eq r2 then // K is not of degree 4; skip printf "ellcdeg4: %o, %o --> smaller degree fields (sqrt(%o), sqrt(%o))\n", pair1, pair2, r1, r2; continue; end if; K := AbsoluteField(ext); X := PolynomialRing(K).1; sqrt1 := Roots(X^2 - rad1)[1,1]; sqrt2 := Roots(X^2 - rad2)[1,1]; parK := [Evaluate(p, X) : p in par]; pol := Denominator(s1)*Denominator(s2)*(parK[2] - sqrt1*parK[3])*(parK[2] - sqrt2*parK[3]) where _, s1 := IsSquare(rad1/r1) where _, s2 := IsSquare(rad2/r2); ptx := Evaluate(xcoord, Eltseq(pt)); polpt := Evaluate(pol, ptx); HE := HyperellipticCurve(polpt*pol); toHE := func; HEcov := map Pr1Q | [HE.1, HE.3]>; Append(~Eseq, ); Append(~info, [pair1, pair2]); end if; end for; // return Eseq; a list of quadruples // with HE a genus 1 hyperelliptic curve over a number field, // toHE maps points on C to sets of points on HE, // HEcov is a morphism from HE to Pr1Q such that relevant points have rational image // (this is needed for the Elliptic Curve Chabauty computation), // xcoord is the rational function on the genus 5 curve giving the x-coordinate on HE return Eseq, info; end function; // Compute an upper bound for the Mordell-Weil rank of an elliptic curve E // over a number field of degree 4. function rank_bound(E, ptsE, Fields) // ptsE is a sequence of points on E; // they provide a lower bound on the rank and can be used to speed up the Selmer group computation // Fields is a set of number fields that will be checked if they can be re-used // (we do not actually use this feature) // // We first compute the lower bound on the rank obtained from the given points printf " computing reduced basis...\n"; time lowbd := #ReducedBasis(ptsE); printf " lower bound from known points = %o\n", lowbd; // Since we are interested in doing Elliptic Curve Chabauty on E, // we need the rank to be less than 4 if lowbd ge 4 then printf " lower bound is too large ==> abort\n"; return 1000, lowbd; end if; // Now comnpute the 2-Selmer group of E to get an upper bound printf " computing 2-Selmer group...\n"; time SelE, mSelE := TwoSelmerGroup(E : Fields := Fields, Hints := {pt[1]/pt[3] : pt in ptsE | pt[3] ne 0}); tors2 := #Invariants(TwoTorsionSubgroup(E)); rkSel := #Invariants(SelE); // upper bound = dim Selmer group - dim E(K)[2] rkbd := rkSel - tors2; printf " rk Sel = %o, rk E[2] = %o ==> rank bound is %o\n", rkSel, tors2, rkbd; if rkbd mod 2 ne lowbd mod 2 then // We do not want to search for missing generators on E. // So if the lower and upper bounds have different parities, then we // stop the computation, since by standard conjectures, the // true rank has the same parity as the upper bound, and so // a gap will remain. printf " parities differ ==> abort\n"; return rkbd, lowbd; end if; if rkbd gt lowbd then // We want the bounds to coincide. If there is a gap, we can look at // 2-isogenous curves (same rank, but possibly better bound). // In our situation, there is always at least one 2-isogenous curve. printf " computing 2-Selmer group of 2-isogenous curve(s)...\n"; T, mT := TwoTorsionSubgroup(E); assert #T ge 2; // can be larger! // each point of order 2 gives a 2-isogenous curve for t in {t : t in T | t ne T!0} do Ei, toEi := IsogenyFromKernel(E, Polynomial([-mT(t)[1], 1])); // tranport known points and saturate at 2 ptsEi := Saturation([toEi(b) : b in ptsE], 2); // compute the 2-Selmer group time SelEi, mSelEi := TwoSelmerGroup(Ei : Hints := {pt[1]/pt[3] : pt in ptsEi | pt[3] ne 0}); tors2 := #Invariants(TwoTorsionSubgroup(Ei)); rkSel := #Invariants(SelEi); // update the upper bound rkbd := Min(rkbd, rkSel- tors2); printf " rk Sel = %o, rk E[2] = %o ==> rank bound is now %o\n", rkSel, tors2, rkbd; end for; end if; return rkbd, lowbd; end function; // Do the Elliptic Curve Chabauty computation. function ECC(E, ptsE, Ecov) // E is the elliptic curve over a number field K, rk E(K) < [K : Q] // ptsE are points on E generating a finite-index subgroup of E(K) // Ecov is a morphism E --> P^1 // This computes the (finite) set of K-points P on E such that Ecov(P) in P^1(Q). bas := ReducedBasis(ptsE); T, mT := TorsionSubgroup(E); lent := #Invariants(T); // we need the (known subgroup of the) Mordell-Weil group of E as an abstract abelian group MW := AbelianGroup(Invariants(T) cat [0 : b in bas]); // and the map from this abstract group to E(K) MWtoE := map E | g :-> mT(T!s[1..lent]) + &+[s[lent+i]*bas[i] : i in [1..#bas]] where s := Eltseq(g)>; // keep track of primes p such that known group is p-saturated // (true for all p up to satbound) satbound := 1; repeat printf " run Chabauty...\n"; // run Elliptic Curve Chabauty time MWset, N := Chabauty(MWtoE, Ecov); // for the result to be provably correct, the known subgroup has to // be saturated at the primes dividing N badp := PrimeDivisors(N); if IsEmpty(badp) or Max(badp) le satbound then // condition is satisfied break; end if; printf " need to saturate at primes in %o\n", [p : p in badp | p gt satbound]; reg := Regulator(bas); // update satbound satbound := Max(badp); // saturate at primes up to (new) satbound time bas := ReducedBasis(Saturation(bas, satbound)); // check if group was enlarged index := Round(Sqrt(reg/Regulator(bas))); if index gt 1 then printf " group enlarged by index %o\n", index; end if; // if group was already saturated, then we are done until index eq 1; printf " group was already saturated\n"; // we return the set of rational points on P^1 that are images of K-points on E xcsE := {Pr1Q!Ecov(MWtoE(g)) : g in MWset}; printf " found %o x-coordinates on E\n", #xcsE; return xcsE; end function; function is_isomorphic(E1, E2); // E1, E2: two elliptic curves over isomorphic base fields that are Galois over Q // check if they are isomorphic over an isomorphism of the base fields; // if so, return also maps E1(K1) -> E2(K2) and E2(K2) -> E1(K1) K1 := BaseField(E1); K2 := BaseField(E2); // get isomorphisms between the base fields flag, K1toK2 := IsIsomorphic(K1, K2); assert flag; autK1, _, mK1 := AutomorphismGroup(K1); auts := [mK1(a) : a in autK1]; j1 := jInvariant(E1); j1inK2 := [K1toK2(a(j1)) : a in auts]; j2 := jInvariant(E2); // for each isomorphism f: K1 --> K2, check if f_* E1 and E2 are isomorphic for i := 1 to #auts do j1K2 := j1inK2[i]; if j1K2 eq j2 then // j-invariants have to be equal E1K2 := BaseChange(E1, auts[i]*K1toK2); flag, iso := IsIsomorphic(E1K2, E2); if flag then return true, func, func; end if; end if; end for; return false, _, _; end function; // This is the main function. // It attempts to find all extensions ("legal" or not) of the given // rational diophantine quadruple aseq = [a_1, a_2, a_3, a_4] to a // rational diophantine quintuple (assuming GRH). function extend(aseq) // some sanity checks assert #aseq eq 4; ChangeUniverse(~aseq, Rationals()); aseq0 := Append(aseq, 0); // "homogenize" the quadruple a,b,c,d := Explode(aseq); // verify that it is really diophantine flag, r12 := IsSquare(a*b+1); assert flag; flag, r13 := IsSquare(a*c+1); assert flag; flag, r14 := IsSquare(a*d+1); assert flag; flag, r23 := IsSquare(b*c+1); assert flag; flag, r24 := IsSquare(b*d+1); assert flag; flag, r34 := IsSquare(c*d+1); assert flag; printf "Trying to find all extensions of the diophantine quadruple\n"; printf " (%o, %o, %o, %o)\n\n", a, b, c, d; // compute the regular extensions abcd := a*b*c*d; part1 := (a+b+c+d)*(abcd+1) + 2*(a*b*c + a*b*d + a*c*d + b*c*d); part2 := 2*r12*r13*r14*r23*r24*r34; ext1 := (part1 + part2)/(abcd-1)^2; ext2 := (part1 - part2)/(abcd-1)^2; assert forall{z : z in aseq | IsSquare(z*ext1 + 1) and IsSquare(z*ext2 + 1)}; printf "regular extensions:\n"; printf " %o and %o\n\n", ext1, ext2; // scale aseq to get coprime integers aseqn := [Integers()!(l*a) : a in aseq] where l := LCM([Denominator(a) : a in aseq]); aseqn := [ExactQuotient(a, g) : a in aseqn] where g := GCD(aseqn); Append(~aseqn, 0); // determine the primes of bad reduction bad := Sort(Setseq(&join[Set(PrimeDivisors(aseqn[j]-aseqn[i])) : j in [i+1..5], i in [1..5]])); printf "bad primes:\n %o\n", bad; // find the five elliptic curves E_i printf "\nconstruct the elliptic curves...\n"; Eseq := ellipticcurves(ChangeUniverse(aseqn, Rationals())); // SelZ is the group of square classes in Q supported on the bad primes SelZ, mSelZ := pSelmerGroup(2, {p*Integers() : p in bad}); // its threefold direct product is the "trivial" upper bound // for the 2-Selmer groups of the E_i SelE, injsE, projsE := DirectProduct([SelZ, SelZ, SelZ]); a,b,c,d := Explode(aseqn); // note roots of right hand sides of the elliptic curves // order is imortant! (must be consistent with iseq and trips below) rtseq := [[-b*c, -b*d, -c*d], [-a*c, -a*d, -c*d], [-a*b, -a*d, -b*d], [-a*b, -a*c, -b*c], [-a*b-c*d, -a*c-b*d, -a*d-b*c]]; assert forall{i : i in [1..5] | Set(rtseq[i]) eq {r[1] : r in Roots(Polynomial([a6,a4,a2,1]))} where a1,a2,a3,a4,a6 := Explode(aInvariants(Eseq[i][2]))}; printf "find their 2-Selmer groups...\n"; // collect information on the groups Sel^(2)(E_i) in imEseq imEseq := []; for i := 1 to 5 do E := Eseq[i,2]; // the curve E_i printf " elliptic curve E_%o:\n", i; // compute 2-Selmer group Sel, mSel := TwoSelmerGroup(E); A := Domain(mSel); // we need to map A to the algebra coming directly from E, // so that we can compute the Selmer group as a subgroup of SelE maxrt := Max(rtseq[i]); scale := (3*maxrt - &+rtseq[i])/Coefficient(Modulus(A), 2); assert scale^3*Evaluate(Modulus(A), (x-maxrt)/scale) eq &*[x - r : r in rtseq[i]]; AtoSelE := map SelE | a :-> &+[injsE[j](mSelZ(Evaluate(f, (rtseq[i,j]-maxrt)/scale))) : j in [1..3]] where f := P!a>; // record the information (the subgroup of SelE corresponding to the Selmer group of E_i) Append(~imEseq, sub); printf " rank bound = %o\n", #Invariants(Sel) - 2; // we could spend more effort here and try to compute the actual group E_i(Q), // but this might take a long time and does not seem to be really necessary end for; printf "\nconstruct image of product of E_i(Q) in H...\n"; // map a point on E to the 2-Selmer group of E function toSelE(pt, rts) if pt[3] eq 0 then return SelE!0; else xs := [pt[1]/pt[3] - r : r in rts]; xs := [xc eq 0 select 0 else sqfp(xc) : xc in xs]; pos := Position(xs, 0); if pos ne 0 then // patch the undefined component so that product is a square xs[pos] := sqfp(&*[xc : xc in xs | xc ne 0]); end if; return &+[injsE[j](mSelZ(xs[j])) : j in [1..3]]; end if; end function; // define the direct sum H_0 of the Sel^(2)(E_i) // (we use this as our initial upper bound for the image of C(Q)) BigH, injsH, projsH := DirectProduct([SelE : i in [1..5]]); H := sub; printf "H_0 has rank %o\n", #Invariants(H); // set up the genus 5 curve given by the quadruple C := genus5curve(aseq); // and the map C(Q) --> H_0 CtoH := map H | pt :-> &+[injsH[i](toSelE(Eseq[i,3](Eseq[i,1]![s[j] : j in [1..5] | j ne i] where s := Eltseq(pt)), rtseq[i])) : i in [1..5]]>; // write down points on C: // points coming from the illegal extension by zero ptsC0 := [C![u1,u2,u3,u4,1] : u1, u2, u3, u4 in [1,-1]]; // points coming from the first regular extension s1, s2, s3, s4 := Explode([s where _, s := IsSquare(1 + aseq[i]*ext1) : i in [1..4]]); ptsC1 := [C![u1*s1,u2*s2,u3*s3,u4*s4,1] : u1, u2, u3, u4 in [1,-1]]; // points coming from the second regular extension s1, s2, s3, s4 := Explode([s where _, s := IsSquare(1 + aseq[i]*ext2) : i in [1..4]]); ptsC2 := [C![u1*s1,u2*s2,u3*s3,u4*s4,1] : u1, u2, u3, u4 in [1,-1]]; // H1 is the image of the trivial points in H_0; it forms a subgroup H1 := {CtoH(pt) : pt in ptsC0}; assert Set(sub) eq H1; // redefine it as a subgroup of H_0 = H H1 := sub; printf "H_1 has rank %o\n", #Invariants(H1); // The action of the automorphism group of C corresponds to translations // by elements of H1, so we can work with the quotient Hq by H1. Hq, toHq := quo; // consider the images of the two "regular" orbits H2a := {CtoH(pt) : pt in ptsC1}; H2b := {CtoH(pt) : pt in ptsC2}; // check H2a, H2b are cosets of H1 assert H2a eq {r + h : h in H1} where r := Rep(H2a); assert H2b eq {r + h : h in H1} where r := Rep(H2b); // note the image of the known points in Hq knownimage := {toHq(Rep(H1)), toHq(Rep(H2a)), toHq(Rep(H2b))}; // look for more rational points on C printf "\nsearching for points on C...\n"; ptsC3 := PointSearch(C, 10^8); ptsC3 := [pt : pt in ptsC3 | pt notin ptset] where ptset := Set(ptsC0) join Set(ptsC1) join Set(ptsC2); if not IsEmpty(ptsC3) then // if extra points were found, add their images to the known image num := ExactQuotient(#ptsC3, 16); printf "found %o additional orbit%o of points on C\n", num, num eq 1 select "" else "s"; for pt in [pt : pt in ptsC3 | forall{i : i in [1..5] | pt[i] ge 0}] do Include(~knownimage, toHq(CtoH(pt))); end for; end if; printf "known points give %o A-orbit%o in H_0\n", #knownimage, #knownimage eq 1 select "" else "s"; // determine a complete set of representatives for the Aut(C)-orbits ptsCrep := []; ptsC := Set(ptsC0 cat ptsC1 cat ptsC2 cat ptsC3); for pt in [pt : pt in ptsC0 cat ptsC1 cat ptsC2 cat ptsC3 | forall{i : i in [1..5] | pt[i] ge 0}] do if pt notin ptsCrep then Append(~ptsCrep, pt); end if; end for; // Compute the hat{\phi}-Selmer set of C: // * produce representatives of the descent functions // * note that scaling aseq scales them by squares, so we can work with aseqn // * also note that we can replace D/N by D*N iseq := [[j : j in [1..5] | j ne i] : i in [1..5]]; Vseq := [&*[aseqn[s[j]] - aseqn[s[i]] : j in [i+1..4], i in [1..4]] : s in iseq]; // Vandermonde P5 := Ambient(C); useq := [u1,u2,u3,u4,u5]; // the D_i (compare Section 3) Dseq := [Determinant(Matrix(Parent(u1), [[1,1,1,1], aseqn[i], [a^2 : a in aseqn[i]], useq[i]])) : i in iseq]; function trips(inds) // coset representatives of Klein Four Group in S_4 return [inds, inds[[1,3,2,4]], inds[[1,4,2,3]]]; end function; // normalize mod squares function norm_sq(pol) cs := ChangeUniverse(Coefficients(pol), Integers()); r, s := SquarefreeFactorization(GCD(cs)); return pol/s^2; end function; // tripseq[i] is the triple of functions on C whose square classes give the map // to the i-th component of BigH tripseq := [[norm_sq(Vseq[i]*Dseq[i]*((useq[t[1]]-useq[t[2]])/(aseqn[t[1]]-aseqn[t[2]]) + (useq[t[3]]-useq[t[4]])/(aseqn[t[3]]-aseqn[t[4]]))) : t in trips(iseq[i])] : i in [1..5]]; // the local version of H at the prime p, with map from BigH function localH(p) assert IsPrime(p); if IsOdd(p) then // Q_p^* mod squares is Z/2Z x Z/2Z Hloc := AbelianGroup([2, 2]); function sqclass(z) v := Valuation(z, p); return Hloc![v, IsSquare(GF(p)!(z/p^v)) select 0 else 1]; end function; else // p = 2: Q_2^* mod squares is Z/2Z x Z/2Z x Z/2Z Hloc := AbelianGroup([2, 2, 2]); function sqclass(z) v := Valuation(z, 2); r := (Integers()!(Integers(8)!(z/2^v))+1) div 2; // 1,3,5,7 --> 1,2,3,4 return Hloc!([v] cat [[0,0],[1,0],[0,1],[1,1]][r]); end function; end if; // construct maps from global to local objects SelZtoloc := hom Hloc | [sqclass(g @@ mSelZ) : g in OrderedGenerators(SelZ)]>; SelEloc, injsEloc, projsEloc := DirectProduct([Hloc, Hloc, Hloc]); SelEtoloc := hom SelEloc | [&+[injsEloc[i](SelZtoloc(projsE[i](g))) : i in [1..3]] : g in OrderedGenerators(SelE)]>; BigHloc, injsHloc, projsHloc := DirectProduct([SelEloc : i in [1..5]]); Htoloc := hom BigHloc | [&+[injsHloc[i](SelEtoloc(projsH[i](g))) : i in [1..5]] : g in OrderedGenerators(H)]>; // "combine" maps a 15-tuple of elements of Hloc to BigHloc combine := func; H1loc := sub; Hqloc, toHqloc := quo; HqtoHqloc := hom Hqloc | [toHqloc(Htoloc(g @@ toHq)) : g in OrderedGenerators(Hq)]>; seqtoHqloc := func; return HqtoHqloc, seqtoHqloc, sqclass; end function; // sign changes coming from A-action signseq := [[1,1,1,1,1], [-1,1,1,1,1], [1,-1,1,1,1], [-1,-1,1,1,1], [1,1,-1,1,1], [-1,1,-1,1,1], [1,-1,-1,1,1], [-1,-1,-1,1,1], [1,1,1,-1,1], [-1,1,1,-1,1], [1,-1,1,-1,1], [-1,-1,1,-1,1], [1,1,-1,-1,1], [-1,1,-1,-1,1], [1,-1,-1,-1,1], [-1,-1,-1,-1,1]]; // obtain local information at the good prime p function local_info_good(p) assert IsPrime(p) and p notin bad; HqtoHqloc, seqtoloc, sqclass := localH(p); result := {Codomain(HqtoHqloc)| }; function evuseq(useq) // find image in Hqloc of a point on C mod p given by its coordinates useq i := 0; repeat // find an orbit representative such that the descent map has good reduction i +:= 1; s := signseq[i]; evseq := [[Evaluate(t, [s[j]*useq[j] : j in [1..5]]) : t in trip] : trip in tripseq]; until forall{trip : trip in evseq | #[t : t in trip | GF(p)!t eq 0] le 1}; for i := 1 to 5 do trip := evseq[i]; tripp := ChangeUniverse(trip, GF(p)); pos := Position(tripp, 0); if pos gt 0 then // fix zero entry using "square norm" condition evseq[i][pos] := &*[t : t in trip | GF(p)!t ne 0]; else assert IsSquare(&* tripp); end if; end for; return seqtoloc(&cat evseq); end function; // points at infinity (i.e. u_5 = 0): they exist iff all a_i are in the same square class if #{sqclass(a) : a in aseqn[1..4]} eq 1 then useq := [Integers()!s where _, s := IsSquare(GF(p)!(aseqn[1]*a)) : a in aseqn]; Include(~result, evuseq(useq)); end if; // other points for z := 0 to p-1 do useq := []; for i := 1 to 4 do flag, s := IsSquare(GF(p)!(1 + aseqn[i]*z)); if not flag then continue z; end if; Append(~useq, Integers()!s); end for; Append(~useq, 1); Include(~result, evuseq(useq)); end for; return HqtoHqloc, result; end function; // keep description of subset of Hq = H/H1 containing image of Selmer set // as union of cosets of a subgroup ker of Hq; // initialize with all of Hq ker := Hq; cosets := {Hq!0}; // include restrictions coming from a new prime function new_ker_cosets(ker, cosets, h, im) // h : Hq --> target, im a subset of target // intersect with h^-1(im) kerh := Kernel(h); newker := ker meet kerh; reps := {z @@ q : z in Q} where Q, q := quo; newcos := {c + r : c in cosets, r in reps | h(c+r) in im}; return newker, newcos; end function; // Now use the good primes up to 100 to cut down the upper bound // on the image of C(Q) in Hq. printf "\nFiltering using good primes up to 100...\n"; for p in PrimesInInterval(3, 100) do if p notin bad then printf " p = %o --> ", p; // get local information Hqtoloc, locinfo := local_info_good(p); // update the upper bound ker, cosets := new_ker_cosets(ker, cosets, Hqtoloc, locinfo); printf "set size = %o * %o = %o\n", #cosets, #ker, #cosets*#ker; if #ker eq 1 and cosets eq knownimage then // reached known lower bound --> stop printf "reduced to image of known points\n"; break; end if; end if; end for; // make sure we know the image (this is not guaranteed, but seems to be very likely) assert #ker eq 1 and cosets eq knownimage; // find representatives for the elements of the Selmer set relpts := [ptsCrep[1]]; rempts := Set(ptsCrep[2..#ptsCrep]); remim := Exclude(cosets, Universe(cosets)!0); while not IsEmpty(remim) do pt := Rep(rempts); Exclude(~rempts, pt); impt := toHq(CtoH(pt)); if impt in remim then Exclude(~remim, impt); Append(~relpts, pt); end if; end while; // determine the sets of points on C in the same class as pt, for each pt in relpts ptsinclass := [{pt1 : pt1 in ptsC | CtoH(pt1) eq impt} where impt := CtoH(pt) : pt in relpts]; // construct F_{T,\xi} printf "\nconstructing the curves F_{T,xi} ...\n"; HEseqs := []; infos := []; // loop over the classes (i.e., the possible \xi's) for i := 1 to #relpts do // construct the curves HEseq, info := ellcdeg4(aseqn, relpts[i]); // add information on the images of the (known) rational points on C HEseq := [*Append(e, &cat[e[2](pt) : pt in ptsC]) : e in HEseq*]; Append(~HEseqs, HEseq); Append(~infos, info); end for; assert #Set(infos) eq 1; // each class should give the same combinatorics infos := infos[1]; // construct the Jacobian elliptic curves of the "genus 1 hyperelliptic curves" F_{T,\xi} printf "constructing their Jacobian elliptic curves...\n"; Eseqs := [[* where E, toE := EllipticCurve(e[1], e[2](relpts[i])[1]) : e in HEseqs[i]*] : i in [1..#relpts]]; if #relpts gt 1 then // check for isomorphisms between the various curves for i := 1 to #infos do for j1 := 1 to #relpts do for j2 := j1+1 to #relpts do // check pairwise E1 := Eseqs[j1,i,1]; E2 := Eseqs[j2,i,1]; flag, isog1, isog2 := is_isomorphic(E1, E2); if flag then printf "found isomorphism: curve %o at points %o and %o\n", i, j1, j2; // this can give us additional known points Eseqs[j2,i,3] cat:= [isog1(pt) : pt in Eseqs[j1,i,3]]; Eseqs[j1,i,3] cat:= [isog2(pt) : pt in Eseqs[j2,i,3]]; end if; end for; end for; end for; end if; // set up a set for collecting the extensions a_5 of (a_1,a_2,a_3,a_4) result := {Rationals()|}; // Now we try to find, for each \xi, one T such that we know // a finite-index subgroup of F_{T,\xi}(K) (as an elliptic curve) rempts := Set(relpts); // classes not yet dealt with successfully // loop over the classes for i := 1 to #relpts do pt := relpts[i]; printf "\ndealing with class of z = %o\n", ((pt[1]/pt[5])^2 - 1)/aseq[1]; // loop over the various elliptic curves F_{T,\xi} for j := 1 to #Eseqs[i] do e := Eseqs[i,j]; he := HEseqs[i,j]; printf " looking at elliptic curve no. %o:", j; for s in infos[j] do printf " %o", s; end for; printf "...\n"; // get rank bounds and hope for equality and < 4 rkbd, lowbd := rank_bound(e[1], e[3], {}); if lowbd lt 4 then if rkbd eq lowbd then // the favorable case; we can do Elliptic Curve Chabauty printf " ==> know finite index subgroup; do ECChab\n"; // transport the map to P^1 from the hyperelliptic to the elliptic model Ecov := Expand(Inverse(e[2])*he[3]); // run Elliptic Curve Chabauty xcsE := ECC(e[1], e[3], Ecov); printf " image on P^1 is %o\n", xcsE; printf " xcoord = %o\n", he[4]; for xcE in xcsE do // check if point lifts to C; xcE is point on Pr1Q CtoP1 := Extend(map Pr1Q | [Numerator(he[4]), Denominator(he[4])]>); lifts := Points(xcE @@ CtoP1); if not IsEmpty(lifts) then ptC := Rep(lifts); z := ((ptC[1]/ptC[5])^2 - 1)/aseq[1]; printf " found extension %o\n", z; Include(~result, z); end if; end for; // now pt is done Exclude(~rempts, pt); continue i; // can skip to next class else printf " rank bound too large\n"; end if; end if; end for; end for; // check if all classes could be dealt with if IsEmpty(rempts) then printf "\nSUCCESS: possible extensions are\n %o\n\n", result; else printf "\nUnable to resolve problem for %o orbits\n\n", #rempts; end if; return result, rempts; end function; // fam(t) gives the quadruple at t in the family considered at the end of Section 3 fam := func; // Now do the computation for certain values of t for t in [2, 3, 2/3, 3/2, 4, 3/4, 4/3, 5, 1/5, 2/5, 3/5, 5/4, 4/5] do printf "\n==================================================================\n"; printf "Considering t = %o\n", t; printf "==================================================================\n\n"; t0 := Cputime(); extensions, remaining := extend(fam(t)); assert IsEmpty(remaining); tt := Cputime(t0); printf "total time: %o s\n", tt; printf "==================================================================\n"; end for;