// 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 [1] 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[1] : e in fdisc | e[2] 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)[2] 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;
assert Dimension(imdelta2) eq g; // sanity check -- independence of the first g basis elements
// Now add the necessary further ones:
d := -1;
while Dimension(imdelta2) lt g + tdim do
d +:= 2;
pol := 4 - x^d;
impol := del2(pol);
if impol notin imdelta2 then
// found a new generator
Append(~basmod2J, pol);
Append(~imbas, impol);
imdelta2 +:= sub
;
end if;
end while;
// Further sanity check: all further x^d - 4 do not give something new.
assert forall {j : j in [d+2..2*g-1 by 2] | del2(4 - x^j) in imdelta2};
// I is the set of d's used for the second bunch of generators.
I := {Degree(pol) : pol in basmod2J[g+1..g+tdim]};
vprintf User1: "I = %o\n\n", I;
// Set up an abstract F_2-vector space to represent imdelta2.
Vimdel2 := KSpace(GF(2), g+tdim);
isoim := hom
;
if Dimension(imU) ne Ngens(U1) then
vprintf User1: "localization map at 2 is not injective\n";
H2group := Codomain(loc);
U1toH2 := hom