// Implement the algorithm for "Selmer Group Chabauty" for curves of the form // y^2 = x^{2g+1} + h(x)^2 // (with some additional conditions on h in Z[x]) // described in "An application of Selmer Group Chabauty' to arithmetic dynamics". // // M. Stoll, started 2019-10-24 // The following is a general function determining the minimum absolute // precision of a logarithm series to achieve a desired 2-adic precision // of the result when evaluating at a suitable value. function necessary_precision(a, b) // a: a positive rational number s.t. the n-th term of the series // has 2-adic valuation >= a*n - v_2(n) + c with some c // b: the desired precision - c // (so that the result is the minimal n such that // (*) a*N - v_2(N) >= b for all N >= n) // First find smallest power of two (>= 1/a) that satisfies (*); // this is an upper bound for n. k := 0; n := 1; while a*n lt Max(1, k+b) do k +:= 1; n *:= 2; end while; // Now do binary search k -:= 1; // 2-adic valuation of the difference of the interval endpoints 2^(k-1), 2^k left := n div 2; right := n; while k gt 0 do middle := (left + right) div 2; k -:= 1; // new interval has half the length; this is also v_2(middle) if a*middle - k ge b then // middle is an upper bound, since a*n - v_2(n) does not get smaller to the right right := middle; else // middle + 1 is a lower bound, since (*) is false for N = middle left := middle; end if; end while; // if k = 0, upper bound is result return right; end function; // The following function determines the number of coefficients of // the series giving the differentials that are necessary to obtain // a given ansolute 2-adic precision for the logarithms, which use // parameters t such that v_2(t) >= 1/g (equality occurring). function necprec(g, prec) // g: genus // prec: target precision // The coefficients a_n in the series omega = SUM_n a_n t^n dt // satisfy v_2(a_n) >= ceil(-2n/(2g+1)). // The series we need to compute have the form // SUM_{n>0} 2^n/n a_{dn-1} (1 <= d <= g) or SUM_{n>0} 4^n/n a_{dn-1} (d = 1,3,...,2g-1). // So we need that // n - v_2(n) + ceil((2-2ng)/(2g+1)) > prec // for n > n0; then result > g n0. // Also, 2n - v_2(n) + ceil((2-2n(2g-1))/(2g+1)) > prec for n > n0'; // then result > (2g-1) n0'. // The first condition is satisfied when // (2+n)/(2g+1) - log_2 n > prec . // When n > (2g+1)/log 2, the lhs is increasing. So start with a reasonable guess // and increase n until the inequality is OK. g2 := 2*g+1; n := g2*(prec + Ceiling(Log(2, g2*prec))); while (2+n)/g2 - Log(2, n) le prec do n +:= 1; end while; result := g*(n-1) + 1; // check second condition n := ((result-1) div (2*g-1)) + 1; assert (2+4*n)/g2 - Log(2, n) gt prec; return result; end function; // This is similar, but for the logarithm series on the Q_2-residue disk around P_0. function necprec1(g, prec) // The condition is // 2/(2g+1) + (2g-1)/(2g+1) n - v_2(n) > prec g2 := 2*g+1; n := Ceiling((g2*prec-2)/(g2-2)); while (2 + (g2-2)*n)/g2 - Log(2, n) le prec do n +:= 1; end while; return n; end function; // The following function checks that the image under \rho log // (i.e. 2 Z_2 -param-> C(Q_2) -i-> J(Q_2) -log*trans-> Z_2^g --> P^2(Q_2) --> P^2(F_2) ) // of a residue disk does not meet what comes from the Selmer group. // The disk is given by a sequence of g power series giving the map // Z_2 -2-> 2 Z_2 --> Z_2^g . function check(logs, prec, T, globinT) // logs: sequence of length g of power series // prec: absolute precision of values of the series when evaluated on Z_2 // T: F_2-vector space of dimension g // globinT: subspace of T that has to be avoided // We do it recursively. // Look at terms of lowest valuation. serprec := Min([AbsolutePrecision(l) : l in logs]) - 1; minval := prec; support := []; for j := 1 to serprec do val := Min([Valuation(Coefficient(l, j), 2) : l in logs]); if val lt minval then minval := val; support := [j]; // new minimum found; elif val eq minval then Append(~support, j); // additional place where minimum is reached end if; end for; error if minval ge prec, "Insufficient precision in formal logs"; // Extract dominant term(s) dom := &+[T!(2^-minval*Vector([Coefficient(l, j) : l in logs])) : j in support]; if dom eq T!0 then // dominant terms cancel --> split into smaller disks vprintf User1: "dom = 0 --> refine\n"; function check_rec(t0, pow2, val, depth) error if val ge prec, "Insufficient precision in formal logs"; dom := T![2^-val*Evaluate(l, t0) : l in logs]; if dom eq T!0 then // recurse vprintf User1: " "^depth*"dom = 0 --> refine\n"; return check_rec(t0, 2*pow2, val+1, depth+1) and check_rec(t0 - Sign(t0)*pow2, 2*pow2, val+1, depth+1); elif dom in globinT then vprintf User1: " "^depth*"dom = %o is in image of global units --> failure!\n", dom; return false; else vprintf User1: " "^depth*"dom = %o is not in image of global units\n", dom; return true; end if; end function; if not check_rec(1, 4, minval+1, 1) or not check_rec(-1, 4, minval+1, 1) then return false; end if; elif dom in globinT then vprintf User1: "dom = %o is in image of global units --> failure!\n", dom; return false; else vprintf User1: "dom = %o is not in image of global units\n", dom; end if; if support eq  then // if linear term is dominant, this won't change for higher valuation of t return true; else // consider next smaller disk return check([Evaluate(l, 2*Universe(logs).1) : l in logs], prec, T, globinT); end if; end function; function SelChab(g, h) // g: natural number >= 1 // h: polynomial in Z[x] h0 := h; // check the arguments error if not (Type(g) eq RngIntElt and g ge 1), "g must be a positive integer"; error if not (Type(h) eq RngUPolElt and Set(Coefficients(h)) subset Integers()), "h must be a polynomial with integral coefficients"; PZ := PolynomialRing(Integers()); h := PZ!h; error if Degree(h) gt g, "The degree of h must be <= g"; error if not (IsOdd(Evaluate(h, 0)) and IsEven(Evaluate(h, 1))), "h(0) must be odd, h(1) must be even"; f := PZ.1^(2*g+1) + h^2; vprintf User1: "SelChab: g = %o, h = %o\n\n", g, h0; disc := Discriminant(f); vprintf User1: "Factoring the discriminant of f ...\n"; vtime User1: fdisc := Factorization(disc); badprimes := {Integers()| e : e in fdisc | e ge 2}; vprintf User1: "S = %o\n", badprimes; // Set up the algebra A P := PolynomialRing(Rationals()); x := P.1; A, toA := quo

; // K as a quotient ring Kprod, AtoKprod := AbsoluteAlgebra(A); // write A as a product of number fields m := NumberOfComponents(Kprod); vprintf User1: "J(Q) has dimension %o\n\n", m-1; vprintf User1: "Computing A(S, 2) ...\n"; vtime User1: U, AtoU, UtoBvec, B := pSelmerGroup(A, 2, badprimes join {2} : Raw); Bseq := Eltseq(B); // sequence of elements in factor base vprintf User1: "Dimension of A(S, 2) is %o\n", #Invariants(U); // Find the subgroup U1 of U compatible with the local conditions at the real place. // Determine the real embeddings and the order of the images of x. r := &+[Signature(Kprod[i]) : i in [1..m]]; // number of real places of A conjth := &cat[Conjugates(imth[i]) : i in [1..m]] where imth := AtoKprod(toA(x)); assert #conjth eq 2*g+1; realplaces := [i : i in [1..2*g+1] | Imaginary(conjth[i]) eq 0]; assert #realplaces eq r; realconjth := [Real(c) : c in conjth[realplaces]]; Sort(~realconjth, ~perm); // sort in ascending order // perm is the permutation such that [Real(conjth[realplaces[i^perm]]) : i in [1..r]] = realconj sortedplaces := [realplaces[i^perm] : i in [1..r]]; prec := Precision(Universe(realconjth)); // Set up H_infty and the map from U into it. Hinf := AbelianGroup([2 : i in [1..r]]); function signs(a) ima := AtoKprod(a); nma := Norm(a); // Compute conjugates until precision is sufficient prec1 := prec; repeat conjs := &cat[Conjugates(ima[i] : Precision := prec1) : i in [1..m]]; prec1 := Ceiling(1.2*prec1); until Abs(nma - &*conjs) lt 1e-10; // now return the signs of the real ones return Hinf![Real(c) gt 0 select 0 else 1 : c in conjs[sortedplaces]]; end function; BinHinf := [signs(b) : b in Bseq]; UtoHinf := hom Hinf | [&+[BinHinf[i] : i in [1..#Bseq] | IsOdd(imu[i])] where imu := UtoBvec(u) : u in OrderedGenerators(U)]>; // The image of J(R) in H_infty (comes from the connected components of C(R)): // subgroup of signs generated by +...+-...- (even number of -'es). allowed := sub; U1 := (allowed meet Image(UtoHinf)) @@ UtoHinf; vprintf User1: "Dimension after local conditions at infinity: %o\n", #Invariants(U1); // Now take into account the local information at the (odd) bad primes. for p in badprimes do locp := LocalTwoSelmerMap(A, p); Hp := Codomain(locp); Qp := pAdicField(p, Max(10, Ceiling(100*Log(2)/Log(p)))); factfp := [e : e in Factorization(PolynomialRing(Qp)!f)]; assert &+[Degree(pol) : pol in factfp] eq 2*g+1; polsp := [P!(-1)^Degree(factfp[i])* (factfp[i] - &*[Universe(factfp)| factfp[j] : j in [1..#factfp] | j ne i]) : i in [1..#factfp]]; imp := sub; if #Invariants(imp) ne #factfp-1 then // Image of 2-torsion is not large enough (<==> nontrivial 4-torsion over Q_p). // Test random potential b-polynomials for points on J over Q_p. Zp := Integers(Qp); count := 1; repeat bpol := Polynomial([Random(Zp) : i in [0..g]]); fb := Factorization(Parent(bpol)!f - bpol^2); for e in fb do if Degree(e) le g then // gives a point on the Jacobian over Q_p imp +:= sub; end if; end for; count +:= 1; until #Invariants(imp) eq #factfp-1 or count gt 1000; if #Invariants(imp) ne #factfp-1 then // Random search unsuccessful --> give up. printf "unable to determine local image at %o\n", p; return false, -p; end if; end if; BinHp := [locp(b) : b in Bseq]; UtoHp := hom Hp | [&+[BinHp[i] : i in [1..#Bseq] | IsOdd(imu[i])] where imu := UtoBvec(u) : u in OrderedGenerators(U)]>; U1 meet:= (imp meet Image(UtoHp)) @@ UtoHp; vprintf User1: "Dimension after local conditions at %o: %o\n", p, #Invariants(U1); end for; // Set up the map delta_2. loc := LocalTwoSelmerMap(A, 2); // this gives the map A^x --> H_2 // The codomain is an elementary abelian 2-group; convert to a vector space. H2 := KSpace(GF(2), #Invariants(Codomain(loc))); // This is the map that one has to apply to (-1)^(deg a) a or variants thereof: del2 := func; // Set up Q_2. Q2 := pAdicField(2, 100); // Determine 2-torsion over Q_2. fact := [e : e in Factorization(ChangeRing(f, Q2))]; tdim := #fact - 1; vprintf User1: "J(Q_2) has dimension %o\n", tdim; // Find image of delta_2 and a-polynomials of generators of J(Q_2)/2J(Q_2). // The image is the image of x^d - 2 (d = 1,...,g) and x^d - 4 (d = 1,3,...,2g-1), // where the first g are independent and we need tdim of the second ones. basmod2J := [(-1)^d*(x^d-2) : d in [1..g]]; imbas := [del2(pol) : pol in basmod2J]; imdelta2 := sub