/******************************************************************** * MAGMA code accompanying the paper * "On the number of rational iterated pre-images of the origin * under quadratic dynamical systems" * by Xander Faber, Benjamin Hutz and Michael Stoll * * September 3, 2010 ********************************************************************/ v1,v2,v3 := GetVersion(); // find out what Magma version is used. if v1 lt 2 or v2 lt 14 then printf "Please use a newer Magma version (2.14 or later).\n"; quit; end if; if v1 eq 2 and v2 lt 16 then printf "This script has been tested on Magma version 2.16;\n"; printf " it also seems to work on version 2.14,\n"; printf " but it may crash version 2.15.\n\n"; end if; // Verification of the proof of Proposition 5.1 printf "\nVerification of the proof of Proposition 5.1\n"; printf "============================================\n\n"; // We work with the projective closure of the affine curve // r^2 = d^2 - d, s^2 = d^2 + d Pr3 := ProjectiveSpace(Rationals(), 3); C := Curve(Pr3, [r^2 - d^2 + d*u, s^2 - d^2 - d*u]); printf "Genus of the curve in Proposition 5.1 is %o\n", Genus(C); // ==> 0, so the curve can be parameterized Pr1 := ProjectiveSpace(Rationals(), 1); printf "Computing parameterization...\n"; para := Parametrization(C, C![1,1,1,0], Curve(Pr1)); F := FunctionField(Rationals()); // extract rational function in t = t1/t0 for d/u (= d in affine chart) dint := Evaluate(DefiningEquations(para)[3], [t,1]) / Evaluate(DefiningEquations(para)[4], [t,1]); printf "Checking that it agrees with the paper (up to sign)...\n"; assert dint eq (t^2+1)^2/(4*t*(t^2-1)); // ==> agrees with paper up to sign cint := -dint^2; // here the sign change is irrelevant preims := [e1*(t^2+1)*(t^2+2*e2*t-1)/(4*t*(t^2-1)) : e1, e2 in [-1,1]]; printf "Second images of x-values: %o\n", [(x^2 + cint)^2 + cint : x in preims]; // ==> [0, 0, 0, 0], shows that we do get four second pre-images // Verification of the proof of Lemma 5.2 printf "\nVerification of the proof of Lemma 5.2\n"; printf "======================================\n\n"; Aff2 := AffineSpace(Rationals(), 2); printf "Defining Y^pre(3,0)...\n"; Ypre30 := Curve(Aff2, ((x^2 + c)^2 + c)^2 + c); eqn30 := DefiningEquation(Ypre30); printf "Its genus is %o\n", Genus(Ypre30); // ==> 1 printf "Defining the elliptic curve E...\n"; E := EllipticCurve([-1,1]); printf "Setting up the maps between E and Y^pre(3,0)...\n"; EtoY := map Ypre30 | [v*q/(u^2-q^2), -q^4/(u^2-q^2)^2]>; // Magma not complaining shows this really defines a map YtoE := map E | [x^2+c, x, (x^2+c)^2+c]>; // dito YtoE := Extend(YtoE); assert IsEmpty(BaseScheme(YtoE)); // ==> this is a morphism EtoY := Extend(EtoY); // this is still undefined at some points EtoE := Extend(Expand(EtoY*YtoE)); // takes a little while printf "The composition E --> Y^pre(3,0) --> E is given by %o;\n", def[#def] where def := AllDefiningPolynomials(EtoE); // ==> [u, v, q], so this composition is the identity printf " since u,v,q are the projective coordinates, this is the identity.\n"; YtoY := Expand(YtoE*EtoY); defpol := DefiningPolynomials(YtoY); // check that this is the identity where it is defined printf "Checking that Y^pre(3,0) --> E --> Y^pre(3,0) is the identity\n"; printf " (where defined)...\n"; assert IsDivisibleBy(Numerator(defpol[1]) - x*Denominator(defpol[1]), eqn30); assert IsDivisibleBy(Numerator(defpol[2]) - c*Denominator(defpol[2]), eqn30); // ==> OK. // look up the curve E printf "The elliptic curve E has Cremona reference %o\n", CremonaReference(E); // ==> 92b1 MWE, toE := MordellWeilGroup(E); printf "Its Mordell-Weil group is\n%o\n", MWE; // ==> free abelian group of rank 1 // find generator genE := toE(MWE.1); printf " and it is generated by P_0 = %o\n", genE; // ==> (1 : 1 : 1), so this is (u, v) = (1, 1). // Verification of the proof of Proposition 5.4 printf "\nVerification of the proof of Proposition 5.4\n"; printf "============================================\n\n"; sik := RealField(30)!SiksekBound(E); // The return value of SiksekBound has the "wrong" type in version 2.14 printf "The Siksek Bound for E is %o\n", ChangePrecision(sik, 5); upper := 4/3*Log(12) + 2*Log(2); lower := -4*sik - 2*Log(2); printf "We get the following bounds:\n"; printf " %o <= 8 h^(P) - h_g(P) <= %o\n", ChangePrecision(lower, 5), ChangePrecision(upper, 5); // note that CanonicalHeight in Magma is with respect to 2*(infty) nbound := Ceiling(((upper-lower)/(4*CanonicalHeight(genE)) - 1)/2); printf "So for n >= %o, the height h_g(n*P_0) is strictly increasing.\n", nbound; // find c-coordinates of points 3*P_0, ..., 20*P_0 printf "Computing c-coodinates of n*P_0, 3 <= n <= 20, ...\n"; clist := [-((pt[1]/pt[3])^2-1)^-2 where pt := n*genE : n in [3..20]]; // and check that their heights are distinct, and the last is the largest hlist := [Max(Abs(Numerator(c)), Denominator(c)) : c in clist]; printf "Checking that they are all distinct and the last has the largest h_g...\n"; assert #Set(hlist) eq #hlist; assert forall{i : i in [1..#hlist-1] | hlist[#hlist] gt hlist[i]}; // Verification of the proof of Proposition 5.5 printf "\nVerification of the proof of Proposition 5.5\n"; printf "============================================\n\n"; if v2 lt 15 then printf "Magma version 2.14: define KillingDifferentialModp...\n"; // Find the reduction mod p of the differential over Q_p that kills // the subgroup generated by a point on the Jacobian function KillingDifferentialModp(pt, p) // Given a point pt of infinite order on the Jacobian of a genus 2 curve // over the rationals of the form y^2 = f(x) and a prime p of good reduction, // this returns a point in P^1(F_p) which is the image of the two zeros // of omega mod p, where omega is a Q_p-differential on the curve that kills // the group generated by pt. J := Parent(pt); C := Curve(J); Jp := BaseChange(J, GF(p)); n := Order(Jp!pt); prec := 1; K := KummerSurface(J); repeat prec +:= 2; K1 := BaseChange(K, pAdicField(p, prec)); ptK := n*K1!K!pt; until Min([Valuation(ptK[i]) : i in [1..3]]) eq prec - 1; coords := [GF(p) | ExactQuotient(Integers()!ptK[i], p^(prec-1)) : i in [1..3]]; assert coords[2]^2 eq 4*coords[1]*coords[3]; Pr1 := ProjectiveSpace(GF(p), 1); return coords[1] eq 0 select Pr1![1,0] else Pr1![coords[2]/2,coords[1]]; end function; end if; printf "Constructing the four curves...\n"; ulist := [x/(x^2+cint) : x in preims]; Clist := [HyperellipticCurve(Numerator((t^2+1)^4*(u^3-u+1))) : u in ulist]; Cpp := Clist[4]; printf "The last is C++:\n%o\n", Cpp; printf "Checking that all curves are isomorphic...\n"; assert forall{C : C in Clist | IsIsomorphic(C, Cpp)}; printf "Defining the genus 2 curve D...\n"; D := HyperellipticCurve(Polynomial([2,2,-1,-3,0,1])); printf "Setting up the morphism C++ --> D...\n"; CpptoD := map D | [-tt^2 - 2*tt*zz + zz^2, (-tt^2 + 2*tt*zz + zz^2)*yy, tt^2 + zz^2]>; // Magma checks that this really defines a map JD := Jacobian(D); rDbound := RankBound(JD); printf "2-descent on the Jacobian of D gives upper bound of %o for the rank\n", rDbound; genJD := D![1,1] - D![-1,-1]; assert Order(genJD) eq 0; // infinite order printf "A point of infinite order is P_D = %o\n", genJD; JDtors, toJDtors := TorsionSubgroup(JD); printf "The torsion subgroup has invariants %o\n", Invariants(JDtors); htc := HeightConstant(JD : Effort := 2); printf "To check that P_D generates the free part of the Mordell-Weil group,\n"; htb := Floor(Exp(CanonicalHeight(genJD)/4 + htc)); printf " we have to verify that all points up to naive height %o\n", htb; printf " are multiples of P_D, up to addition of a torsion point.\n"; ptsJD := Points(JD : Bound := htb); ptsJD1 := {n*genJD + toJDtors(t) : n in [-6..6], t in JDtors}; assert forall{pt : pt in ptsJD | pt in ptsJD1}; printf "This is the case.\n"; printf "Now do the Chabauty argument.\n"; printf "Check that 3 is a good prime...\n"; assert 3 notin BadPrimes(D); kdiff := KillingDifferentialModp(genJD, 3); printf "The reduction of the differential over Q_3 killing the MW group\n"; printf " vanishes mod 3 at the points with x-coordinate %o\n", kdiff; printf "These points are not defined over F_3.\n"; assert not IsSquare(Evaluate(PolynomialRing(GF(3))!HyperellipticPolynomials(D), kdiff[1]/kdiff[2])); printf "This implies that there is at most one rational point on D\n"; printf " in each residue class mod 3.\n"; numptsD3 := #BaseChange(D, Bang(Rationals(), GF(3))); printf "There are %o residue classes,\n", numptsD3; ptsD := Points(D : Bound := 1000); printf " and we know the following %o points on D:\n%o\n", #ptsD, ptsD; assert numptsD3 eq #ptsD; printf "So these must be all the rational points on D.\n"; printf "Lifting these back to C++, we obtain the set of rational points on C++:\n"; ptsCpp := &join[Points(pt @@ CpptoD) : pt in ptsD]; assert forall{pt : pt in ptsCpp | pt[1]*pt[3]*(pt[1]^2-pt[3]^2) eq 0}; printf "%o\nThere are %o points, all with t in {0, -1, 1, infty}.\n", ptsCpp, #ptsCpp; if v1 gt 2 or v2 ge 15 then printf "Note:\n"; printf "The built-in Magma function Chabauty(JacHypPt) does a combined\n"; printf " Chabauty and Mordell-Weil sieve computation and determines D(Q):\n"; printf " D(Q) = %o\n", Chabauty(genJD); else printf "Note:\n"; printf "As of version 2.15,\n"; printf " the built-in Magma function Chabauty(JacHypPt) does a combined\n"; printf " Chabauty and Mordell-Weil sieve computation and determines D(Q).\n"; end if; // some helper functions // test a point for regularity on an arithmetic surface function isregular(S, pt) // S: affine or projective curve over Z // pt: point in S(F_p) IS := Ideal(S); P := Generic(IS); r := Rank(P); p := #Universe(Eltseq(pt)); Iloc := ideal

; if IsProjective(S) then // pick suitable affine chart j := r; while pt[j] ne 1 do j -:= 1; end while; assert j gt 0; IS := IS + ideal

; end if; I2 := IS + Iloc^2; dim := Dimension(ChangeRing(quo

, GF(p))); return dim eq 2 // point is regular on fiber or (p in I2 and dim eq 3); // if p in I2 then OK, else dimension is too large end function; // find primes of bad reduction function BadPrimesArithmeticSurface(S) // S: projective curve over Z SS := SingularSubscheme(S); IS := Ideal(SS); P := Generic(IS); r := Rank(P); ISseq := [ideal

: Z in [P.j : j in [1..r]]]; // affine patches ISelimseq := [EliminationIdeal(I, {}) : I in ISseq]; // intersect with Z badfact := Factorization(LCM([Integers()!Basis(I)[1] : I in ISelimseq])); return Sort([a[1] : a in badfact]); end function; // Verification of the proof of Theorem 6.2 printf "\nVerification of the proof of Theorem 6.2\n"; printf "========================================\n\n"; printf "Defining X^pre(4,0) in P^4 over the integers...\n"; PZ5 := PolynomialRing(Integers(), 5); Pr4Z := ProjectiveSpace(PZ5); XZ := Curve(Pr4Z, [Z3^2+Z1*Z4-Z0^2, Z3^2+Z2*Z4-Z1^2, Z3^2+Z3*Z4-Z2^2]); bad := BadPrimesArithmeticSurface(XZ); printf "Its bad primes are %o\n", bad; eqns := DefiningEquations(XZ); PEF := PolynomialRing(Integers()); printf "Computing its Euler factor over F_3...\n"; EF3 := Numerator(ZetaFunction(BaseChange(XZ, GF(3)))); printf "Its factorization is\n%o\n", Factorization(EF3); printf "Computing the Picard groups of J^pre(4,0) over F_p\n"; printf " for p = 3, 5, 7, 11 and 29...\n"; plist := [3, 5, 7, 11, 29]; assert IsEmpty(Set(plist) meet Set(bad)); results := []; for p in plist do Gp, mGp := ClassGroup(BaseChange(XZ, GF(p))); Append(~results, [* Gp, mGp *]); printf " p = %o done\n", p; end for; torsbd := 0; for pair in results do torsbd := GCD(torsbd, &*[Integers() | i : i in Invariants(pair[1]) | i ne 0]); end for; printf "Torsion bound (GCD of group orders) is %o\n", torsbd; printf "Since Pic^0 over F_3 is cyclic: %o\n", [i : i in Invariants(results[1][1]) | i ne 0]; printf " and Pic^0 over F_29 is killed by 2*odd: %o\n", [i : i in Invariants(results[5][1]) | i ne 0]; printf " it follows that the torsion subgroup has order at most 2.\n"; XQ := BaseChange(XZ, Rationals()); ptsXQ := [XQ | [1,1,1,1,0], [1,1,1,-1,0], [1,1,-1,1,0], [1,1,-1,-1,0], [1,-1,1,1,0], [1,-1,1,-1,0], [1,-1,-1,1,0], [1,-1,-1,-1,0], [0,0,0,0,1], [0,-1,0,-1,1] ]; printf "The list of known rational points on X^pre(4,0):\n%o\n", ptsXQ; places := [Place(pt) : pt in ptsXQ]; relations := [[-1, 0, 1, 1, 0,-1, 1, 1,-1,-1], [-1, 1, 1,-1,-1, 1, 1,-1, 0, 0], [-1,-1, 1, 0,-1,-1, 0, 1, 1, 1], [-2,-1, 0, 0, 1, 1, 1, 1,-1, 0], [-1,-1, 0, 1, 1, 0,-1,-1, 2, 0], [ 1, 1,-1, 0, 0,-1, 1, 1, 0,-2]]; printf "Checking that the following divisors are principal:\n"; for r in relations do divisor := &+[r[i]*places[i] : i in [1..10]]; sf := true; printf " "; for j := 1 to 10 do if r[j] ne 0 then if r[j] gt 1 then printf sf select "%o*P_%o" else " + %o*P_%o", r[j], j; elif r[j] eq 1 then printf sf select "P_%o" else " + P_%o", j; elif r[j] lt -1 then printf sf select "-%o*P_%o" else " - %o*P_%o", -r[j], j; else // r[j] eq -1 printf sf select "-P_%o" else " - P_%o", j; end if; sf := false; end if; end for; printf " ==> %o\n", IsPrincipal(divisor); end for; printf "This shows that the rank of the subgroup G generated by\n"; printf " degree-zero divisors supported in P_1, ..., P_10\n"; printf " has rank at most (10-1) - 6 = 3.\n"; printf "Setting up map from Pic/Q to product of Pic/F_p...\n"; reductions := [ [* pair[1], [(1*Place(Xp!ChangeUniverse(Eltseq(pt), GF(p)))) @@ pair[2] : pt in ptsXQ] *] where Xp := Curve(Codomain(pair[2])) where pair := results[i] where p := plist[i] : i in [1..#plist] ]; FA10 := FreeAbelianGroup(10); Picprod, incls, projs := DirectProduct([pair[1] : pair in reductions]); FA10toPicprod := hom Picprod | [&+[incls[i](reductions[i][2][j]) : i in [1..#plist]] : j in [1..10]]>; degzero := Kernel(hom Z | [Z|[1] : i in [1..10]]>) where Z := FreeAbelianGroup(1); imdegzero := FA10toPicprod(degzero); printf "Invariants of image of degree zero part of known group are\n%o\n", Invariants(imdegzero); printf "Since this group has (Z/7Z)^3 as a quotient,\n"; printf " this implies that G has rank 3.\n"; Pic0, qPic0 := quo; printf "Checking that generates G...\n"; assert Image(hom Pic0 | [qPic0(FA10.2-FA10.3), qPic0(FA10.2-FA10.5), qPic0(FA10.2-FA10.8)]>) eq Pic0; // Verification of the proof of Theorem 6.3 printf "\nVerification of the proof of Theorem 6.3\n"; printf "========================================\n\n"; p := 3; printf "Chabauty: working at p = %o\n", p; Xp := BaseChange(XZ, GF(p)); ptsp := Points(Xp); printf "Found %o points mod %o\n", #ptsp, p; printf "Checking that known rational points surject onto points mod %o...\n", p; imptsQ := {@ Xp!ChangeUniverse(Eltseq(pt), GF(p)) : pt in ptsXQ @}; assert imptsQ eq ptsp; printf "Verify kernel of reduction at %o...\n", p; basindex := 5; basept := ptsXQ[basindex]; printf "Choosing %o as base-point\n", basept; FA10toClp := hom reductions[1][1] | reductions[1][2]>; Pic0toClp := hom reductions[1][1] | [FA10toClp(g @@ qPic0) : g in OrderedGenerators(Pic0)]>; Rkerreps := [FA10 | [0,1,0,1,-6,1,1,0,1,1], [0,-1,0,-1,2,-1,1,-1,0,1], [2,-1,0,0,-3,2,0,-1,0,1]]; Rkergens := [qPic0(r) : r in Rkerreps]; assert Kernel(Pic0toClp) eq sub; printf "Points Q_1, Q_2, Q_3 as in paper generate intersection of MW group\n"; printf " with kernel of reduction at %o\n", p; printf "Changing to a plane model (this is less bug-prone)...\n"; X1, toX1 := Projection(XQ, XQ![0,0,0,0,1]); X2, toX2 := Projection(X1, X1![0,1,0,1]); X2 := Curve(X2); printf "Plane model:\n%o\n", X2; printf "Constructing the map to the new model...\n"; time prmap := Extend(Restriction(Expand(toX1*toX2), XQ, X2)); printf "Matching up the places on the old and new models...\n"; ptsX2 := {@ prmap(pt) : pt in ptsXQ @}; assocpl := []; f := FunctionField(XQ)!((P.1-P.2)/P.2) where P := CoordinateRing(Ambient(XQ)); ff := Pushforward(prmap, f); // use f to tell points apart for ptX2 in ptsX2 do ptsabove := [ : i in [1..#ptsXQ] | prmap(ptsXQ[i]) eq ptX2]; if #ptsabove eq 1 then Append(~assocpl, ); else plc := Places(ptX2); assert #ptsabove eq 2 and #plc eq 2; assert Valuation(f, ptsabove[1,2]) ne Valuation(f, ptsabove[2,2]); if Valuation(ff, plc[1]) eq Valuation(f, ptsabove[1,2]) then assocpl cat:= [, ]; else assocpl cat:= [, ]; end if; end if; end for; Sort(~assocpl, func); assert [a[2] : a in assocpl] eq ptsXQ; places2 := [a[3] : a in assocpl]; printf "Places on new model (same order as points on old model):\n%o\n", places2; printf "Finding representative divisors for kernel of reduction...\n"; basDiv := 1*places2[basindex]; repDivs := [Reduction(&+[s[j]*places2[j] : j in [1..#s]] where s := Eltseq(r), basDiv) : r in Rkerreps]; X2p := BaseChange(X2, Bang(Rationals(), GF(p))); printf "Taking affine patch with last coordinate = 1\n"; X2A := AffinePatch(X2, 1); P2Q := CoordinateRing(Ambient(X2A)); h := DefiningEquation(X2A); printf "Affine patch given by\n%o\n", h; bpt := X2A!RepresentativePoint(places2[basindex]); printf "Base-point on this affine patch is %o\n", bpt; assert Valuation(Evaluate(Derivative(h, 2), Eltseq(bpt)), p) eq 0; pwsprec := 50; Pws := LaurentSeriesAlgebra(Rationals(), pwsprec); uint := t + bpt[1]; uniformizer := u - bpt[1]; printf "Can take t = %o as uniformizer at base-point\n", uniformizer; printf "Solve for v as power series in t = %o ...\n", uniformizer; solns := Roots(Evaluate(h, [uint, P.1])) where P := PolynomialRing(Pws); solns := [r[1] : r in solns | Coefficient(r[1], 0) eq bpt[2]]; assert #solns eq 1; vint := solns[1]; printf " ==> v = %o\n", vint + BigO(t^6); printf "Extracting representative points of the divisors...\n"; decomps := [Decomposition(D) : D in repDivs]; assert forall{dec : dec in decomps | #dec eq 1 and dec[1,2] eq 1}; divplcs := [dec[1,1] : dec in decomps]; printf "Finding the characteristic polynomial of their t-values...\n"; charpolys := [CharacteristicPolynomial(RepresentativePoint(pl)[1] - bpt[1]) : pl in divplcs]; assert forall{j : j in [1..#divplcs] | Degree(charpolys[j]) eq Degree(divplcs[j])}; printf "Checking that points all reduce to the base-point mod %o...\n", p; assert forall{pol : pol in charpolys | forall{j : j in [0..Degree(pol)-1] | Valuation(Coefficient(pol, j), p) gt 0}}; printf "Convert differentials on X^pre(4,0) into power series in t, times dt...\n"; Pr4Q := Ambient(XQ); P5Q := CoordinateRing(Pr4Q); FF2Q := FieldOfFractions(P2Q); // pull back uniformizer tonXQ := Pullback(prmap, FunctionField(X2A)!uniformizer); dtonXQ := Differential(tonXQ); DB0 := DifferentialBasis(Divisor(XQ, Scheme(Pr4Q, P5Q.1))); assert #DB0 eq 1; omega0 := DB0[1]; omegafns := [FF2Q!Pushforward(prmap, omega0/dtonXQ)]; for j := 2 to Rank(P5Q) do Append(~omegafns, omegafns[1]*FF2Q!Pushforward(prmap, FunctionField(XQ)!(P5Q.j/P5Q.1))); end for; omegas := [Evaluate(f, [uint, vint]) : f in omegafns]; printf "Found basis corresponding to hyperplanes z0 = 0, ..., z4 = 0:\n"; for j := 1 to #omegas do printf " omega_%o = (%o) dt\n", j-1, omegas[j] + O(t^5); end for; logs := [Integral(om) : om in omegas]; printf "The corresponding logarithms (= integrals) are\n"; for j := 1 to #logs do printf " l_%o = %o\n", j-1, logs[j] + O(t^6); end for; printf "Computing the power sums of the t-values for the generators...\n"; powersums := [[Coefficient(s, i) : i in [1..pwsprec]] where s := -t*Derivative(f)/f where f := &+[Coefficient(pol, Degree(pol)-j)*t^j : j in [0..Degree(pol)]] : pol in charpolys]; // now do the integration intprec := 5; Qpprec := 2*intprec; cofreq := 5*(intprec - 1); e := 0; repeat e +:= 1; n := p^e*Floor(5*(intprec+e-1)/p^e); cofreq := Max(cofreq, n); until n eq 0; assert cofreq lt Min([AbsolutePrecision(l) : l in logs]); Qp := pAdicField(p, Qpprec); printf "Evaluating the integrals to precision O(%o^%o)...\n", p, intprec; mat := Matrix([[Qp!&+[Coefficient(l, i)*s[i] : i in [1..cofreq]] + O(Qp!p^intprec) : s in powersums] : l in logs])/p; printf " ==> this results in (%o times) the matrix\n%o\n", p, Transpose(mat); matker1 := Basis(Kernel(ChangeRing(mat, pAdicRing(p, intprec-1)))); matker := Basis(Kernel(ChangeRing(mat, GF(p)))); assert #matker eq 2 and #matker1 eq 2; printf "Its kernel is spanned (mod %o^%o) by\n%o\n", p, intprec-1, matker1; printf "The kernel of its reduction mod %o is spanned by\n%o\n", p, matker; P5p := CoordinateRing(Ambient(Xp)); hyperplanes := [&+[r[j]*P5p.j : j in [1..Rank(P5p)]] : r in matker]; printf "This corresponds to the pencil spanned by the hyperplanes\n"; printf " %o = 0 and %o = 0 in P^4/F_%o\n", hyperplanes[1], hyperplanes[2], p; printf "Now checking that for every point P_j, j /= 1,8, in X^pre(4,0)(F_%o),\n", p; printf " there is a differential in the pencil that does not vanish at P\n"; badpoints := [Parent()| ]; for pt in ptsp do printf " point %o ==> ", pt; if exists{h : h in hyperplanes | Evaluate(h, Eltseq(pt)) ne 0} then printf "OK\n"; else j := 1; while pt[j] eq 0 do j +:= 1; end while; val := Min([Valuation(FunctionField(Xp)!(h/P5p.j), pt) : h in hyperplanes]); printf "both differentials vanish: minimal valuation is %o\n", val; Append(~badpoints, ); end if; end for; PwsQp := LaurentSeriesAlgebra(Qp, pwsprec); // find global points corresponding to the bad points above badindex := [i : i in [1..#ptsXQ] | Xp!ChangeUniverse(Eltseq(ptsXQ[i]), GF(p)) in [a[1] : a in badpoints]]; // check that the bad points are the reductions of P_1 and P_8 assert badindex eq [1, 8]; badplaces := places2[badindex]; AFF, toAFF := AlgorithmicFunctionField(FunctionField(X2A)); PF := PolynomialRing(GF(p)); for pl in badplaces do pt := X2A!RepresentativePoint(pl); printf "\nLooking at `bad' place %o...\n", pl; plFF := FunctionFieldPlace(pl); // take completion printf "Computing completion of function field at this place...\n"; Loc, toLoc := Completion(AFF, plFF); printf "Expressing differentials in terms of the uniformizer at this place...\n"; uniformizeratpt := UniformizingParameter(pl); factoratpt := toAFF((Differential(FunctionField(X2A)!uniformizer) / Differential(uniformizeratpt))); omegasatpt := [toLoc(toAFF(f)*factoratpt) : f in omegafns]; for j := 1 to #omegasatpt do printf " omega_%o = (%o) dt\n", j-1, omegasatpt[j] + O(t^5); end for; logsatpt := [Integral(om) : om in omegasatpt]; printf "The corresponding logarithms are\n"; for j := 1 to #logsatpt do printf " l_%o = %o\n", j-1, logsatpt[j] + O(t^6); end for; killinglogsatpt := [&+[b[j]*PwsQp!logsatpt[j] : j in [1..5]] : b in matker1]; killinglogsatpt := [Evaluate(l, p*T) : l in killinglogsatpt]; valmins := [ where m, pos := Min([Valuation(c) : c in Coefficients(l)]) : l in killinglogsatpt]; minprecs := [Min([AbsolutePrecision(c) : c in Coefficients(l)]) : l in killinglogsatpt]; cutoffs := [Max([j : j in [1..AbsolutePrecision(l)-1] | Valuation(Coefficient(l, j)) lt minprecs[i]]) where l := killinglogsatpt[i] : i in [1..#minprecs]]; rich := [ : i in [1..#cutoffs]]; klatpt := [&+[(Coefficient(a[1],j) + O(Qp!p^a[3]))/Coefficient(a[1],a[2,2]) * T^j : j in [1..a[4]]] + O(T^(a[4]+1)) : a in rich]; printf "The functions on the residue class vanishing on rational points are\n"; for ser in klatpt do printf " %o\n", ser + O(T^intprec); end for; printf "Extracting mod %o parts gives\n", p; polys := [PF!ChangeUniverse([Coefficient(l, j) : j in [0..AbsolutePrecision(l)-1]], GF(p)) : l in klatpt]; for pol in polys do printf " %o\n", pol; end for; polgcd := GCD(polys); printf "Their gcd is %o\n", polgcd; assert Degree(polgcd) eq 1; printf " degree is 1 ==> only one rational point in residue class\n"; end for; // Verification of the local information in Section 7 printf "\nVerification of the local information in Section 7\n"; printf "==================================================\n\n"; _ := PEF; assert bad eq [ 2, 23, 2551 ]; p := bad[3]; // 2551 printf "Looking at reduction mod %o\n", p; Xp := BaseChange(XZ, GF(p)); SXp := SingularSubscheme(Xp); assert Degree(SXp) eq 1; singptp := Rep(Points(SXp)); printf "Unique singular point at %o\n", singptp; printf "Checking regularity...\n"; assert isregular(XZ, singptp); printf " ==> given model is regular at %o\n", p; assert #Places(Xp!singptp) eq 2; printf "Tangents at singularity split\n"; numptsp1 := #Points(Xp) + 1; // add one since have two points on the smooth curve // EulerFactor is (1 - T)(1 - (p+1-numptsp1) T + ...) EF2551 := 1 -(1 + (p+1-numptsp1))*T; printf "Euler factor at %o starts %o + ...\n", p, EF2551; f2551 := 1; printf "Conductor exponent at %o is %o\n", p, f2551; printf "\n"; p := bad[2]; // 23 printf "Looking at reduction mod %o\n", p; Xp := BaseChange(XZ, GF(p)); SXp := SingularSubscheme(Xp); deg := Degree(SXp); singptsp := Points(SXp); assert #singptsp eq deg and deg eq 3; printf "Three singularities at\n"; for j := 1 to deg do printf " %o", singptsp[j]; if j lt deg then printf ","; else printf "\n"; end if; end for; printf "Checking regularity...\n"; for pt in singptsp do assert isregular(XZ, pt); end for; printf " ==> given model is regular at %o\n", p; for sing in singptsp do assert #Places(Xp!sing) eq 2; end for; printf "Tangents at all singularities split\n"; assert IsIrreducible(Xp); // ==> genus 2 curve with three nodes // Euler factor is (1 - T)^3 times Euler factor of smooth genus 2 curve EF23 := (1 - T)^3*PEF!Numerator(ZetaFunction(Xp)); printf "Euler factor at %o is %o\n", p, EF23; f23 := 3; printf "Conductor exponent at %o is %o\n", p, f23; printf "\n"; p := bad[1]; // 2 printf "Looking at reduction mod %o\n", p; Xp := BaseChange(XZ, GF(p)); SXp := SingularSubscheme(Xp); deg := Degree(SXp); assert deg eq 16; assert Degree(ReducedSubscheme(SXp)) eq 1; singptp := Xp!Rep(Points(SXp)); plc := Places(singptp); assert #plc eq 1 and Degree(plc[1]) eq 1; printf "Unique singular point at %o (a higher cusp);\n", singptp; assert not isregular(XZ, singptp); printf " this point is not regular on the arithmetic surface.\n"; printf "Moving to an affine model in A^3 with singularity at origin\n"; printf " (similar to the model in the proof of Lemma 7.4)\n"; P4Z := PolynomialRing(Integers(), 4); eqns1 := [Evaluate(e, [z0, z1, z2, 1, z2^2-1]) : e in eqns]; assert eqns1[3] eq 0; eqns1 := [Evaluate(eqns1[i], [z0+z1+z2+1, z1+z2+1, z2+1, z3]) : i in [1..2]] cat [z3-2]; eqns1[1] -:= eqns1[2]; printf "Equations of new model:\n"; for e in eqns1 do printf " %o = 0\n", e; end for; printf "We blow up the singular point:\n"; printf " z0 <-- y0*w, z1 <-- y1*w, z2 <-- y2*w, z3 <-- y3*w\n"; printf " where (y0:y1:y2:y3) are projective coordinates,\n"; printf " and we divide the equations by the highest possible power of w.\n"; P5Z := PolynomialRing(Integers(), 5); printf "First, we replace 2*z0, 2*z1 by z0*z3, z1*z3:\n"; eqns2 := [eqns1[1] - z0*eqns1[3], eqns1[2] - z1*eqns1[3], eqns1[3]]; for e in eqns2 do printf " %o = 0\n", e; end for; eqnsbl := [ExactQuotient(Evaluate(eqns2[1], [w*y0,w*y1,w*y2,w*y3]), w^2), ExactQuotient(Evaluate(eqns2[2], [w*y0,w*y1,w*y2,w*y3]), w^2), Evaluate(eqns2[3], [w*y0,w*y1,w*y2,w*y3])]; printf "This gives the new equations\n"; for e in eqnsbl do printf " %o = 0\n", e; end for; polmod2 := func; printf "The reduction mod 2 is given by\n"; for e in eqnsbl do printf " %o = 0\n", polmod2(e); end for; printf "The component with y3 = 0 is the (strict transform of the) original curve.\n"; printf "The new components are those with w = 0.\n"; printf "We see four lines (in P^3):\n"; printf " y0 = y1 = 0, y0 = y1 + y3 = 0, y0 + y3 = y1 = 0, y0 + y3 = y1 + y3 = 0\n"; printf "They all intersect at (w,y0,y1,y2,y3) = (0,0,0,1,0).\n"; printf "We dehomogenize: (w,y0,y1,y2,y3) <-- (z0,z1,z2,1,z3). This gives\n"; eqnsbl1 := [Evaluate(e, [z0,z1,z2,1,z3]) : e in eqnsbl]; for e in eqnsbl1 do printf " %o = 0\n", e; end for; assert not isregular(Scheme(AffineSpace(P4Z), eqnsbl1), [GF(2)|0,0,0,0]); printf "The origin is not regular.\n"; printf "We replace 2, 2*z1, 2*z2 by z0*z3, z0*z1*z3, z0*z2*z3:\n"; eqnsbl2 := [eqnsbl1[1] + (z2-z1)*eqnsbl1[3], eqnsbl1[2] + eqnsbl1[3], eqnsbl1[3]]; for e in eqnsbl2 do printf " %o = 0\n", e; end for; printf " and blow up again. This gives\n"; eqnsblnew := [ExactQuotient(Evaluate(eqnsbl2[1], [w*y0,w*y1,w*y2,w*y3]), w^2), ExactQuotient(Evaluate(eqnsbl2[2], [w*y0,w*y1,w*y2,w*y3]), w), Evaluate(eqnsbl2[3], [w*y0,w*y1,w*y2,w*y3])]; for e in eqnsblnew do printf " %o = 0\n", e; end for; printf "The reduction mod 2 is given by\n"; for e in eqnsblnew do printf " %o = 0\n", polmod2(e); end for; printf "y3 = 0 gives the original curve:\n"; for e in eqnsblnew[[1,2]] do printf " %o = 0\n", Evaluate(polmod2(e), [w,y0,y1,y2,0]); end for; printf "y0 = 0 (and w /= 0 generically) gives the four lines from the last blow-up:\n"; for e in eqnsblnew[[1,2]] do printf " %o = 0\n", Evaluate(polmod2(e), [w,0,y1,y2,y3]); end for; printf "whereas w = 0 gives the new components:\n"; for e in eqnsblnew[[1,2]] do printf " %o = 0\n", Evaluate(polmod2(e), [0,y0,y1,y2,y3]); end for; printf "This gives the two lines\n"; printf " w = y0 = y1 = 0 and w = y0 = y1+y3 = 0,\n"; printf "each of multiplicity 3 (2 from w^2 = 0, 1 from y0 = 0).\n"; printf "They meet each other and the original curve at\n"; printf " (w,y0,y1,y2,y3) = (0,0,0,1,0)\n"; printf "The first meets two of the other lines at (0,0,0,0,1) and (0,0,0,1,1);\n"; printf " the second meets the other two lines at (0,0,1,0,1) and (0,0,1,1,1).\n"; printf "We check that all the points on the two mutiple lines are regular.\n"; printf "Modulo ^2, the equations are\n"; function redmodsq(pol, vars) ms := Monomials(pol); cs := Coefficients(pol); return &+[(cs[i] mod 2)*ms[i] : i in [1..#ms] | &+[Degree(ms[i], v) : v in vars] le 1]; end function; for e in eqnsblnew do printf " %o = 0\n", redmodsq(e, [w,y0,y1]); end for; printf "The corresponding matrix is\n"; matjac := Matrix([[polmod2(Evaluate(Derivative(e, v), [0,0,0,y2,y3])) : v in [w,y0,y1]] : e in eqnsblnew[[1,2]]]); printf "%o\n", matjac; minors := Minors(matjac, 2); printf "The relevant minors are\n"; for e in minors do printf " %o\n", polmod2(e); end for; assert GCD([Evaluate(m, [0,0,0,PP.1,PP.2]) : m in minors] where PP := PolynomialRing(GF(2), 2)) eq 1; printf "They do not all vanish on P^1/F_2.\n"; printf "This shows that all points are regular,\n"; printf " since the matrix always has rank 2.\n"; printf "For the other line, we first substitute y1 <-- y1-y3.\n"; printf "Then the line is again given by w = y0 = y1 = 0,\n"; printf " and, proceeding as above, we obtain the minors\n"; matjac := Matrix([[polmod2(Evaluate(Derivative(e, v), [0,0,0,y2,y3])) : v in [w,y0,y1]] where e := Evaluate(ee, [w,y0,y1-y3,y2,y3]) : ee in eqnsblnew[[1,2]]]); minors := Minors(matjac, 2); for e in minors do printf " %o\n", polmod2(e); end for; assert GCD([Evaluate(m, [0,0,0,PP.1,PP.2]) : m in minors] where PP := PolynomialRing(GF(2), 2)) eq 1; printf "As before, they do not all vanish simultaneously.\n"; // Verification of the proof of Lemma 7.3 printf "\nVerification of the proof of Lemma 7.3\n"; printf "======================================\n\n"; printf "We substitute (Z0,Z1,Z2,Z3,Z4) <-- (1,z1+1,z2+1,z3+1,z4)\n"; printf " and include the equation z5^3 = 2 (standing for z5 = 2^(1/3))\n"; printf "This gives\n"; PolZ5 := PolynomialRing(Integers(), 5); eqnst3a := [Evaluate(e, [1,z1+1,z2+1,z3+1,z4]) : e in eqns]; eqnst3a := [eqnst3a[1], eqnst3a[2]-eqnst3a[1], eqnst3a[3]-eqnst3a[1]] cat [z5^3 - 2]; for e in eqnst3a do printf " %o = 0\n", e; end for; printf "Now we replace 2 by z5^3 in the first three equations:\n"; eqnst3b := [eqnst3a[1] + z3*eqnst3a[4], eqnst3a[2] - z1*eqnst3a[4], eqnst3a[3] - z2*eqnst3a[4], eqnst3a[4]]; for e in eqnst3b do printf " %o = 0\n", e; end for; printf "The next step is to substitute\n"; printf " z1 <-- (x1*x3+x2)^2 * x2^4 * x3^3 * x4^3\n"; printf " z2 <-- (x1*x3+x2)^2 * x2^3 * x3^3 * x4^2\n"; printf " z3 <-- (x1*x3+x2) * x2^2 * x3^2 * x4^2\n"; printf " z4 <-- (x1*x3+x2)^2 * x2^4 * x3^4 * x4^3 * x5\n"; printf " z5 <-- (x1*x3+x2) * x2 * x3 * x4\n"; printf "and to divide by suitable powers of (x1*x3+x2), x2, x3, x4.\n"; PolZ5x := PolynomialRing(Integers(), 5); subs := [(x1*x3+x2)^2 * x2^4 * x3^3 * x4^3, (x1*x3+x2)^2 * x2^3 * x3^3 * x4^2, (x1*x3+x2) * x2^2 * x3^2 * x4^2, (x1*x3+x2)^2 * x2^4 * x3^4 * x4^3 * x5, (x1*x3+x2) * x2 * x3 * x4]; eqnst3c := [Evaluate(e, subs) : e in eqnst3b]; eqnst3c := [ExactQuotient(eqnst3c[1], (x1*x3+x2)^2 * x2^4 * x3^4 * x4^3), ExactQuotient(eqnst3c[2], (x1*x3+x2)^4 * x2^7 * x3^6 * x4^5), ExactQuotient(eqnst3c[3], (x1*x3+x2)^3 * x2^6 * x3^6 * x4^4), eqnst3c[4]]; eqnst3c[2] := ExactQuotient(eqnst3c[2] - x3*eqnst3c[1] - x2*x4*eqnst3c[4], x3*x4); printf "This gives the model\n"; for e in eqnst3c do printf " %o = 0\n", e; end for; AffF2 := AffineSpace(GF(2), 5); special := Scheme(AffF2, eqnst3c); irr := IrreducibleComponents(special); printf "Its reduction mod 2 has %o irreducible components:\n", #irr; for comp in irr do flag := IsReduced(comp); comp1 := comp; if not flag then comp1 := ReducedSubscheme(comp); end if; for b in MinimalBasis(comp1) do printf " %o =", b; end for; printf " 0;\n (geometric) genus = %o\n", Genus(Curve(comp1)); end for; // Verification of the proof of Lemma 7.4 printf "\nVerification of the proof of Lemma 7.4\n"; printf "======================================\n\n"; printf "Setting up Z[pi], pi = 2^(1/7), ...\n"; R := ext; pi := R!R.2; assert pi^7 eq 2; PR := PolynomialRing(R, 3); eqnst7 := [Evaluate(e, [z0,z1,z2,2]) : e in eqns1[[1,2]]]; PR1 := PolynomialRing(GF(2), 3); red := hom PR1 | hom PR1 | 0>, [x0,x1,x2]>; printf "After substituting\n"; printf " (z0,z1,z2) <-- (pi^7*x0, pi^6*x1, pi^4*x2)\n"; printf " and dividing by suitable powers of pi, we obtain mod pi:\n"; eqnst7a := [Evaluate(e, [pi^7*z0, pi^6*z1, pi^4*z2]) : e in eqnst7]; exquo := func; eqnst7a := [exquo(eqnst7a[1], pi^14), exquo(eqnst7a[2], pi^12)]; for e in eqnst7a do printf " %o = 0\n", red(e); end for; special := Curve(AffineSpace(PR1), [red(e) : e in eqnst7a]); assert IsIrreducible(special) and IsReduced(special) and Genus(special) eq 3; printf "This defines a curve of (geometric) genus 3.\n"; EF2 := PEF!1; printf "\nEuler factor at %o is %o\n", 2, EF2; f2 := 10; printf "Conductor exponent at %o is %o\n", 2, f2; // Computing the L-series values printf "\nComputation of the values of the L-series and its derivatives\n"; printf "=============================================================\n\n"; cond := 2^f2*23^f23*2551^f2551; printf "Conductor %o =", cond; for pf in Factorization(cond) do printf " %o", pf[1]; if pf[2] gt 1 then printf "^%o", pf[2]; end if; end for; printf " ...\n"; prec := 4; // precision bound for L-series computation, higher value takes longer Ls := LSeries(2, [0,0,0,0,0,1,1,1,1,1], cond, 0 : Sign := -1, Precision := prec); // expect odd rank, so sign = -1 numcofs := LCfRequired(Ls); printf "For a precision of %o decimal digits, need %o L-series coefficients\n", prec, numcofs; // For a later check that the sign is correct: Ls1 := LSeries(2, [0,0,0,0,0,1,1,1,1,1], cond, 0 : Sign := 1, Precision := prec); // primes to deal with for the coefficients pr4 := [p : p in [3..Floor(numcofs^(1/4)) by 2] | IsPrime(p) and not p in bad]; pr3 := [p : p in [pr4[#pr4]+2..Floor(numcofs^(1/3)) by 2] | IsPrime(p) and not p in bad]; pr2 := [p : p in [pr3[#pr3]+2..Floor(numcofs^(1/2)) by 2] | IsPrime(p) and not p in bad]; pr1 := [p : p in [pr2[#pr2]+2..numcofs by 2] | IsPrime(p) and not p in bad]; printf "Need to compute Euler factors\n"; printf " to degree %o or higher for %o <= p <= %o\n", 4, Min(pr4), Max(pr4); for j := 3 to 1 by -1 do prl := [pr1,pr2,pr3,pr4][j]; printf " to degree %o for %o <= p <= %o\n", j, Min(prl), Max(prl); end for; // use a more direct approach to point counting on X^pre(4,0) function countpoints40(F) // count F-rational points on X^pre(4,0), F a finite field of char /= 2, 23, 2551 count := 8; // points at infinity for d in F do // (d,-d^2), (-d,-d^2) are preimages of 0; only need one, since -d is also used // note: if d /= -d, then preimages of d will be distinct from those of -d c := -d^2; set := {F| }; flag1, s := IsSquare(d - c); if flag1 then // s, -s are second preimages flag2, t := IsSquare(s - c); if flag2 then // t, -t are third preimages flag3, u := IsSquare(t - c); if flag3 then set join:= {u, -u}; end if; flag3, u := IsSquare(-t - c); if flag3 then set join:= {u, -u}; end if; end if; flag2, t := IsSquare(-s - c); if flag2 then // t, -t are third preimages flag3, u := IsSquare(t - c); if flag3 then set join:= {u, -u}; end if; flag3, u := IsSquare(-t - c); if flag3 then set join:= {u, -u}; end if; end if; end if; count +:= #set; end for; return count; end function; EFs1 := [<2551, EF2551>]; EFs2 := []; EFs3 := []; EFs4 := [<2, EF2>, <23, EF23>]; printf "Computing Euler factors for high degree (using ZetaFunction)...\n"; for p in pr4 do Append(~EFs4, ); printf "p = %2o done\n", p; end for; function coeffs3(p, l) // converts l = [N_1,N_2,N_3] where N_j = #C(\F_p^j) into a,b,c such that // Euler factor starts 1 + a T + b T^2 + c T^3 a := l[1] - (p+1); b := ExactQuotient(l[2] - (p^2+1) + a^2, 2); c := ExactQuotient(l[3] - (p^3+1) - a^3 + 3*a*b, 3); return a, b, c; end function; printf "Computing Euler factors for degree 3...\n"; for p in pr3 do counts := [countpoints40(GF(p,n)) : n in [1..3]]; Append(~EFs3, ); printf "p = %3o done\n", p; end for; printf "Computing Euler factors for degree 2...\n"; for p in pr2 do counts := [countpoints40(GF(p,n)) : n in [1..2]]; Append(~EFs2, ); printf "p = %3o done\n", p; end for; printf "Computing Euler factors for degree 1...\n"; lastp := 1; for p in pr1 do Append(~EFs1, ); if Floor(p/1000) gt Floor(lastp/1000) then printf "%6o done\n", p; end if; lastp := p; end for; EFs := EFs4 cat EFs3 cat EFs2 cat EFs1; EFmap := pmap PEF | EFs>; // Set coefficient function LSetCoefficients(Ls, func); LSetCoefficients(Ls1, func); printf "Checking functional equation...\n"; time f1 := CheckFunctionalEquation(Ls); time f2 := CheckFunctionalEquation(Ls1); printf " ==> %o for sign -1, %o for sign +1\n", f1, f2; printf "This indicates that our choice of sign is correct\n"; printf "(To get a really convincingly small value,\n"; printf " one has to use higher precision...)\n"; // Now compute L(1), L'(1), L''(1), L'''(1) printf "Now we evaluate the L-series and its derivatives at s = 1.\n"; printf " evaluating L(X^pre(4,0), 1)...\n"; time val0 := Evaluate(Ls, 1); printf " ==> %o\n", val0; printf " evaluating L'(X^pre(4,0), 1)...\n"; time val1 := Evaluate(Ls, 1 : Derivative := 1, Leading); printf " ==> %o\n", val1; printf " evaluating L''(X^pre(4,0), 1)...\n"; time val2 := Evaluate(Ls, 1 : Derivative := 2, Leading); printf " ==> %o\n", val2; printf " evaluating L'''(X^pre(4,0), 1)...\n"; time val3 := Evaluate(Ls, 1 : Derivative := 3, Leading); printf " ==> %o\n", val3; printf " ==> rank is bounded above by 3 (assuming BSD+)\n\n";