// Determine the rational points on y^2 = x^6 + x^2 + 1
// ====================================================
//
// This comes from a problem from Diophantus (the only one ivolving
// a higher genus curve):
// ==> Find (all) rational numbers a such that the square of a,
// the square of the square of a and the square of that add up to a square,
// i.e.: a^2 + a^4 + a^8 = b^2. Division by a^2 leads to the equation above.
//
// Joe Wetherell in his thesis solved the problem, using Chabauty's
// method on covering curves.
// We will use Chabauty's method on elliptic curves instead,
// a method pioneered by Nils Bruin.
//
// Michael Stoll, Goettingen, July 2006
//
// NOTE: The code below works with version 2.13-1
//
//======================================================================
//
// Set up polynomial ring
P := PolynomialRing(Rationals());
// Set up curve
f := x^6 + x^2 + 1;
C := HyperellipticCurve(f);
// Look for points
ptsC := Points(C : Bound := 1000); ptsC;
// It is likely that this is all. The only interesting solutions are
// x = +-1/2, y = +-9/8.
//
// Let us first try to see if we can use Chabauty directly.
J := Jacobian(C);
time RankBounds(J);
// So the rank is 2, which is just too large.
// There are two elliptic curves C maps to:
// y^2 = x^3 + x + 1 and y^2 = x^3 + x^2 + 1
// If one of them has rank zero, we are also fine.
RankBounds(EllipticCurve([1,1]));
RankBounds(EllipticCurve([0,1,0,0,1]));
// Bad luck...
//
// We can factor f over K = Q(w) where w^3 + w + 1 = 0:
K := NumberField(x^3+x+1);
PK := PolynomialRing(K);
Factorization(PK!f);
g := $1[1,1]; h := $1[2,1];
// If P = (xx,yy) is a rational point on C, then g(xx) = d u^2 for some
// d in K(S, 2), where S is the set of prime ideals of K dividing the
// resultant of g and h.
OK := Integers(K);
Factorization(ideal);
// This computes K(S, 2):
H, mH := pSelmerGroup(2, {$1[1,1]}); H;
// We want to find the image of C(Q) in H.
// It is contained in the image of J(Q) under the map that on points on C
// is given by (xx,yy) |--> xx^2 - w.
// We need to find (an odd index subgroup of) J(Q).
// What we get from the known points on C is the following.
gens := ReducedBasis([J![ptsC[i], ptsC[1]] : i in [2..#ptsC]]); gens;
// The image of the first point in H is trivial,
// the image of the second point is -w.
// We have to check that our subgroup is of odd index.
// First verify there is no torsion:
TorsionSubgroup(J);
// One way of checking our subgroup is 2-saturated is to look at reductions
// mod suitable primes p. For example:
J3 := BaseChange(J, GF(3));
AbelianGroup(J3);
Order(J3!gens[1]), Order(J3!gens[2]), Order(J3!(gens[1]+gens[2]));
// This shows that at most gens[1] can be divisible by 2.
// In fact, the image of gens[1] in J(F_p) will always be divisible by 2,
// for any prime p of good reduction.
// So we have to use a different method: we try to divide gens[1] by 2.
// This is best done using the Kummer Surface.
Kum := KummerSurface(J);
// Turn it into a scheme:
Pr3<[z]> := ProjectiveSpace(Rationals(), 3);
KumS := Scheme(Pr3, DefiningPolynomial(Kum));
// Set up the duplication map:
dup := map KumS | ChangeUniverse(Kum`Delta, CoordinateRing(Pr3))>;
// Find the image of gens[1] on the Kummer Surface:
pt1 := KumS!Eltseq(Kum!gens[1]);
// Preimage scheme under dup:
pt1 @@ dup;
// Any rational points in there?
Points($1);
// No, so gens[1] is not divisible by 2 in J(Q),
// and our group is 2-saturated.
// (Note: if there were points, one still would have to check if they pull
// back to rational points on J...)
//
// We can conclude that the only relevant elements of H are 0 and -w.
mH(-w);
// So -w is nontrivial in H, and every rational point on C gives rise
// to a K-rational point on y^2 = d h(x) with rational x-coordinate,
// where d = 1 or d = -w.
// These are genus 1 curves, and since we have a point on each of them,
// they are in fact elliptic curves.
E1a := HyperellipticCurve(h);
E2a := HyperellipticCurve(-w*h);
E1, mE1 := EllipticCurve(E1a); E1;
// This does not work:
E2, mE2 := EllipticCurve(E2a);
// but we know a point (with x-coordinate zero):
Points(E2a, 0);
E2, mE2 := EllipticCurve(E2a, $1[1]); E2;
// Now we need to find the MW groups of these elliptic curves.
time success, A1, A1toE1 := PseudoMordellWeilGroup(E1); success; A1;
time success, A2, A2toE2 := PseudoMordellWeilGroup(E2); success; A2;
// So the second curve has rank zero and only two points.
// These must be the two points on E2a with x-coordinate zero,
// hence they come from the points (0, +-1) on C.
//
// For E1, we have to use Chabauty (possible since rank = 1 < 3 = [K:Q]).
// We need the x-coordinate map E1a -> P^1 on E1.
Pr1 := ProjectiveSpace(Rationals(), 1);
pr := map Pr1 | [R.1, R.3] where R := CoordinateRing(Ambient(E1a))>;
prE1 := Extend(Expand(Inverse(mE1)*pr));
// Now we can apply Chabauty:
time n, v, r := Chabauty(A1toE1, prE1, 5 : Aux := {7});
n;
// There are at most 6 K-points on E1
// that on E1a have rational x-coordinate.
v;
// These are the elements of the MW group of E1 giving these points.
r;
// The statements above are valid if index of the given group in MW is
// prime to r. Since r = 2 and PseudoMordellWeilGroup produces something
// of odd index, we are fine.
//
// Let us find the relevant points on C.
ptsE1 := {A1toE1(pt) : pt in v}; ptsE1;
// Magma still has problems pulling back points to weighted P^2's,
// so we take the points we know and check whether this is what we get.
ptsE1a := &join{Points(E1a, xx) : xx in {Infinity(), 1/2, -1/2}}; ptsE1a;
mE1 := Extend(mE1); // to make it defined everywhere
ptsE1 eq {mE1(pt) : pt in ptsE1a};
//
// We have now proved:
//
// THEOREM. The only rational solutions to y^2 = x^6 + x^2 + 1
// are (x, y) = (0, +-1), (+-1/2, +-9/8).