// Arithmetic of Hyperelliptic Curves in MAGMA
// Michael Stoll
// ICTP, September 8, 2017
//
// Define the polynomial ring in the variable x over the rationals.
P := PolynomialRing(Rationals());
//
// Set up one of the two curves of genus 2 with the largest number
// of known rational points in the LMFDB (as of Wednesday).
C := HyperellipticCurve(-2*x^4 + x^3 + 7*x^2 + 4*x, x^3 + x + 1); C;
//
// We can enumerate its rational points up to a given height bound,
// but first we have to change to a "simplified model" y^2 = f(x).
C1, CtoC1 := SimplifiedModel(C);
ptsC1 := Points(C1 : Bound := 10^4); ptsC1; #ptsC1;
// We can move the points back to C:
ptsC := {@ pt @@ CtoC1 : pt in ptsC1 @};
//
// The discriminant of C is prime:
disc := Discriminant(C); disc;
Factorisation(Integers()!disc);
//
// This implies that its conductor is prime as well (the same prime):
Conductor(C);
//
// Let us try to determine the Mordell-Weil group.
// Set up the Jacobian (we use C1):
J1 := Jacobian(C1);
//
// We first check that the torsion subgroup is trivial:
TorsionSubgroup(J1);
//
// The following does a 2-Selmer group computation
// (and possibly a bit more) to bound the rank:
RankBound(J1);
//
// So let us see if we can find that many independent points.
// We can use the points we found on the curve:
pt0 := ptsC1[1]; // base-point for the embedding into the Jacobian
ptsJ1 := [pt - pt0 : pt in ptsC1]; ptsJ1;
bas, mat := ReducedBasis(ptsJ1); bas; mat;
// So the rank is 4, and J1(Q) = Z^4.
//
// One possibility to saturate the known subgroup is to find all
// points of canonical height up to the square of the covering radius
// of the corresponding lattice (with h^ as norm).
// We can get an upper bound for this square as 1/4 times
// the sum of the norms of the basis of the orthogonalised lattice:
ortho := OrthogonalizeGram(mat);
crbound := &+[ortho[i,i] : i in [1..4]]/4; crbound;
//
// To do the enumeration, we have to bound the naive height.
// We get a bound for the difference between naive and canonical height
// as follows:
diffbd := HeightConstant(J1 : Effort := 2); diffbd;
//
// So we have to search for points up to (logarithmic) naive height
bound := crbound + diffbd; bound;
// corresponding to multiplicative naive height
hbound := Floor(Exp(bound)); hbound;
//
// Do the search:
ptsJ66 := Points(J1 : Bound := hbound); #ptsJ66;
// restrict to the points whose canonical height is below the bound:
time {pt : pt in ptsJ66 | Height(pt) le crbound};
// check that we already know the points:
assert $1 subset {J1!0} join &join{{b,-b} : b in bas};
//
// So the points in bas are already generators of J1(Q).
//
// Preview: The following will eventually be included in Magma.
// You can download the code from
// http://www.mathe2.uni-bayreuth.de/stoll/magma/index.html
// under "Computation of the Mordell-Weil group of Jacobians of genus 2";
// then the following four lines of code can be carried out.
//
// AttachSpec("~/magma/g2wrapper/g2wrapper.spec");
// time MordellWeilGroupGenus2(J1 : Quiet := false);
// [m(g) : g in OrderedGenerators(G)] where G, m := $1;
// bas;
//
// Let us verify BSD (numerically) for J1 (or J = Jac(C)).
//
// Set up the L-series (precision not too large for time reasons):
Ls := LSeries(C : Precision := 10);
LCfRequired(Ls); // how many coefficients do we need?
//
// Check the order of vanishing at s=1
time Evaluate(Ls, 1.0);
time Evaluate(Ls, 1.0 : Derivative := 1, Leading);
time Evaluate(Ls, 1.0 : Derivative := 2, Leading);
time Evaluate(Ls, 1.0 : Derivative := 3, Leading);
time Evaluate(Ls, 1.0 : Derivative := 4, Leading);
// the first four values are zero to our working precision,
// whereas the last one clearly is nonzero --> OK!
//
bsd := $1;
// Conjecturally,
// bsd/4! = #Sha(J) Reg(J(Q)) Omega(J) C(J)/#J(Q)_tors^2
// The torsion subgroup is trivial.
// Since v_p(disc) = 1 for the only prime p dividing the discriminant,
// the product of Tamagawa numbers C(J) is 1.
// #Sha should be a square:
IsEven(J1);
// (this also follows from the fact that C(Q) is non-empty).
//
// We compute the regulator:
reg := Determinant(mat); reg;
//
// To get the real period Omega(J), we use the analytic Jacobian.
Jan := AnalyticJacobian(C1);
// The information is contained in its big period matrix.
bpm := BigPeriodMatrix(Jan);
// We need the area of a fundamental domain for the lattice in R^2
// generated by the traces (from C to R) of the columns of bpm.
realbpm := Matrix([[2*Real(bpm[i,j]) : j in [1..4]] : i in [1..2]]);
Matrix(2, 4, [ChangePrecision(z, 10) : z in Eltseq(realbpm)]);
// The last two columns are zero, so
omega := Abs(Determinant(Submatrix(realbpm, [1..2], [1..2]))); omega;
//
// (Side note: since C is regular over Z, the relevant basis of
// differentials is the standard one,
// dx/(2y + x^3+x+1) and x dx/(2y + x^3+x+1).
// The big period matrix is computed w.r.t. dx/y and x dx/y on C1.
// Fortunately, the difference in scaling is exactly cancelled by
// the scaling induced by changing the model from C to C1.)
//
// Now check:
sha_an := bsd/Factorial(4)/reg/omega; sha_an;
// OK!
// (We know at least that Sha[2] = 0, so #Sha should be an odd square.)
//
//=====================================================================
//
// Using Chabauty's method + Mordell-Weil Sieve
//
// Set up a curve (not in the LMFDB)
C := HyperellipticCurve(4*x^6 - 12*x^5 - 23*x^4 - 54*x^3 + 217*x^2 + 408*x + 144);
// It actually has a better model at p=2:
Cred := ReducedMinimalWeierstrassModel(C); Cred;
Factorisation(Integers()!Discriminant(Cred));
Factorisation(Conductor(Cred));
//
// Its Jacobian has rank 1:
J := Jacobian(C);
RankBound(J);
//
// We find 16 rational points on C (record for rank 1):
ptsC := Points(C : Bound := 10^4); ptsC; #ptsC;
//
// We now check that these are all.
// We need a point of infinite order in J(Q).
// The Jacobian has large rational torsion:
Invariants(TorsionSubgroup(J));
// so we may have to search a bit.
[Order(pt - ptsC[1]) : pt in ptsC];
ptJ := ptsC[Position($1, 0)] - ptsC[1]; ptJ; Order(ptJ);
//
// Now we call the Chabauty function.
chab, primes, run := Chabauty(ptJ); chab;
chab eq Set(ptsC);
//
// How 'Chabauty' works:
// * It finds a good odd prime p such that the p-adic differential
// that kills J(Q) under the integration pairing has reduction mod p
// and does not vanish at any F_p-point of C.
// * This implies that there is at most one rational point
// in each residue class mod p on C.
// * Then it uses the "Mordell-Weil Sieve" to show that there are
// no rational points in some residue classes.
// * At the same time, it searches for points in the residue classes.
// * As soon as for each residue class, it has either found a point
// or proved that the class does not contain rational points, it is done.