// Determine rational n-cycles under z |--> z^2 + c
// ================================================
//
// I.e., pairs of rational numbers (z,c) such that
// with z_0 = z, z_{m+1} = z_m^2 + c,
// we have z_0, z_1, ..., z_{n-1} distinct and z_n = z_0.
//
// For example, (0, 0) gives a rational 1-cycle,
// (0, -1) gives a rational 2-cycle 0 |-> -1 |-> 0, ...
//
// See Flynn, Poonen and Schaefer: Cycles of quadratic polynomials ...
// Duke Math. J. 90 (1997).
//
// Michael Stoll, Göttingen, Tue, December 14, 2004
//
// NOTE: The code below works with version 2.11-11
//
// Set up polynomial ring in two variables z,c
P2zc := PolynomialRing(Rationals(), 2);
//
// Find the first few iterates of f : z |--> z^2 + c
iterates := [i eq 0 select z else Self(i)^2 + c : i in [0..6]];
iterates[1];
iterates[2];
iterates[3]; // and so on.
//
// The condition that z_n = z_0 = z means iterates[n+1] - z = 0.
// So for 1-cycles, we get the relation between z and c given by
g1 := iterates[2] - z; g1;
//
// This obviously describes a conic in the affine (z,c)-plane,
// which we can parametrize by z |--> (z, z-z^2)
//
// For 2-cycles, we have iterates[3] = z, but this is also satisfied
// for 1-cycles, so we have to divide by g1.
g2 := ExactQuotient(iterates[3] - z, g1); g2;
//
// This is again a conic, which is trivially parametrized by
// z |--> (z, -z^2-z-1).
//
// For 3-cycles, we get in the same way
g3 := ExactQuotient(iterates[4] - z, g1); g3;
//
// This is not so easily recognized, so we use Magma's plane curves machinery.
// We first have to set up the ambient space, here it is an affine plane.
Aff2 := AffineSpace(P2zc);
//
// Then we turn our equation into an affine curve.
C3 := Curve(Aff2, g3); C3;
//
// Now we can ask Magma for its genus.
Genus(C3);
//
// So despite its complex appearance, it is a genus 0 curve.
// If it has a rational point, then we can again parametrize all the
// rational points.
// In order to do that, it is best to first transform it into a conic.
// Before we can do that, however, we need a projective model.
C3p := ProjectiveClosure(C3);
Pr2 := Ambient(C3p);
C3p;
// Turn this into a conic:
Can := CanonicalDivisor(C3p);
phi := DivisorMap(-Can);
C3i := Image(phi); C3i;
phi := Restriction(phi, C3p, C3i); // to have correct codomain
// We need its inverse:
_, phii := IsInvertible(phi);
Co3, phi1 := Conic(C3i);
// Check if the conic has a rational point:
HasPoint(Co3);
// Parametrize it:
L__ := ProjectiveSpace(Rationals(), 1);
phi2 := Parametrization(L, Co3);
// Compose maps to get parametrization of C3:
// par := Expand(phi2*phi1^-1*phii); // problem here (fixed by Nils Bruin)
phi3 := Expand(phi1^-1*phii);
def := [Evaluate(f, pols) : f in DefiningPolynomials(phi3)]
where pols := DefiningPolynomials(phi2);
def;
// Make them a bit nicer:
def := [ExactQuotient(f, g) : f in def] where g := GCD(def);
def := [f/g : f in def]
where g := GCD(ChangeUniverse(&cat[Coefficients(f) : f in def],
Integers()));
def;
// With some more work, you get it down to this:
def :=
[-2*u*v*(u - v)*(u^3 - u*v^2 + v^3),
-u^6 + 2*u^5*v - 4*u^4*v^2 + 8*u^3*v^3 - 9*u^2*v^4 + 4*u*v^5 - v^6,
4*u^2*v^2*(u - v)^2];
par := map C3p | def>;
// Note that the fact that Magma accepts this input shows that the
// polynomials in def really define a parametrization of C3p!
// Now we can generate examples...
[par(L![a,b]) : a, b in [-2..2] | GCD(a,b) eq 1];
{C3!pt : pt in $1 | pt[3] ne 0};
//
// Let us move on to n=4. Here, we have to remove 1-cycles and 2-cycles
// from the equation.
g4 := ExactQuotient(iterates[5] - z, g1*g2); g4;
//
// What is its genus?
C4 := Curve(Aff2, g4);
time Genus(C4);
//
// So this is a genus 2 curve and therefore must have a model of the form
// y^2 = F(x) with F of degree 5 or 6.
//
// How can we find it?
// We first need a basis of the regular differentials on C4.
Om4, mOm4 := SpaceOfHolomorphicDifferentials(C4);
omega1 := mOm4(Om4.1); // first differential
omega2 := mOm4(Om4.2); // second differential
//
// For every such basis, there is a model y^2 = F(x) such that
// omega1 = dx/y and omega2 = x dx/y.
// So x = omega2/omega1 and y = dx/omega1.
fx := omega2/omega1;
fy := Differential(fx)/omega1;
fx;
// fy fills a screen...
//
// We have to find the relation between fx and fy.
// One way of doing that is to use them to define a map into Pr2 and
// then look at its image.
// But first, we have to convert fx and fy into quotients of polynomials.
fxn := P2zc!Numerator(fx, C4); fxd := P2zc!Denominator(fx, C4);
fyn := P2zc!Numerator(fy, C4); fyd := P2zc!Denominator(fy, C4);
// (here it turns out that fx and fy are already polynomials...)
//
// Now set up the map.
mC4 := map Pr2 | [fxn/fxd, fyn/fyd, 1]>;
time immC4 := Image(mC4);
immC4;
//
// This gives us a nice hyperelliptic curve.
// Now we have to define it as a hyperelliptic curve in Magma.
P1 := PolynomialRing(Rationals());
H4 := HyperellipticCurve(Evaluate(DefiningPolynomial(immC4), [x,0,1]));
H4;
// (This happens to be X_1(16)...)
//
// Now set up the map from C4 to H4.
// NOTE: H4 lives in a weighted projective plane with weights (1,3,1).
// (Here it doesn't matter since den=1).
// Also note that in the equation we got above, we had 9y^2 and not y^2.
mC4 := map H4 | [fxn, 3*fyn, 1]>;
//
// To answer the original question, we need to find the set of rational
// points on H4.
//
// We can first do a simple search for points.
ptsH4 := Points(H4 : Bound := 1000); ptsH4;
//
// What are the corresponding points on C4?
[ pt @@ mC4 : pt in ptsH4 ]; // ==> error message related to weighted P^2
// So we have to do this in a more pedestrian way, by solving algebraic
// equations.
[Variety(ideal) : pt in ptsH4 | pt[3] ne 0];
//
// so there are no _affine_ points on C4 that map to the points we found
// on H4.
//
// It remains to prove that the list of points on H4 is already complete.
//
// To do this, we will use information on the Mordell-Weil group of
// the Jacobian of H4. So we first construct the Jacobian.
J4 := Jacobian(H4);
//
// Next, we have to find the rank of the MW group.
// A bound for the rank is obtained through 2-descent.
SetVerbose("JacHypSelmer", 2); // to see what's going on...
time TwoSelmerGroupData(J4);
// The first value is zero; this proves that the MW rank is zero.
// (In general, this gives an upper bound on the rank.)
//
// So J4(Q) is just the torsion subgroup. Let us determine what it is.
TorsionSubgroup(J4);
//
// Now, using a suitable map of H4 into J4, we can easily find the set
// of rational points on H4. It is still simpler to use a function
// provided for that purpose.
Chabauty0(J4);
$1 eq ptsH4;
// So we have indeed found all points and proved the following.
// THEOREM. There are no rational 4-cycles under z |--> z^2 + c.
//
// What about n = 5?
//
g5 := ExactQuotient(iterates[6] - z, g1); g5;
C5 := Curve(Aff2, g5);
time Genus(C5);
// This took a while...
//
// This genus is too large to allow direct computations with C5.
//
// So what can we do?
//
// Note that we have an action of the cyclic group Z/5Z on C5, given by
// the map we iterate: (z, c) |--> (z^2+c, c).
// We can hope to get a curve we can more easily deal with by dividing
// by this action.
// We find the quotient by setting up a suitable map and asking for its
// image.
mC5 := map Aff2 | [&+[iterates[i] : i in [1..5]], c]>;
D5 := Image(mC5); D5;
Genus(D5);
// So we are back in business!
//
// We repeat the procedure we did for C4 to find a hyperelliptic model.
Om5, mOm5 := SpaceOfHolomorphicDifferentials(D5);
omega1 := mOm5(Om5.1); // first differential
omega2 := mOm5(Om5.2); // second differential
//
fx := omega2/omega1;
fy := Differential(fx)/omega1;
//
fxn := P2zc!Numerator(fx, D5); fxd := P2zc!Denominator(fx, D5);
fyn := P2zc!Numerator(fy, D5); fyd := P2zc!Denominator(fy, D5);
//
mD5 := map Pr2 | [fxn/fxd, fyn/fyd, 1]>;
time immD5 := Image(mD5);
immD5;
//
// Define a hyperelliptic curve; be careful to take the coefficient
// of y^2 into account.
def := DefiningPolynomial(immD5);
H5 := HyperellipticCurve(-1/MonomialCoefficient(def, Y^2*Z^4)
*Evaluate(def, [x,0,1]));
H5;
//
// and the map to it.
mD5 := map H5 | [fxn/fxd, fyn/fyd, 1]>;
//
// Now the model H5 is not particularly nice,
// so let us find a better one.
M5, mH5 := ReducedMinimalWeierstrassModel(H5); M5;
//
// For further compuations, it is better to remove the y terms.
S5, mM5 := SimplifiedModel(M5); S5;
//
// Now we can again look for points.
ptsS5 := Points(S5 : Bound := 1000); ptsS5;
//
// Let us try to move them back to D5.
ptsM5 := {@ pt @@ mM5 : pt in ptsS5 @}; ptsM5;
ptsH5 := {@ pt @@ mH5 : pt in ptsM5 @}; ptsH5;
//
// Now we would run into the same bug as before, so use the work-around:
// (note that all points are affine, with pt[3] = 1)
h5 := DefiningPolynomial(D5);
[Variety(ideal) : pt in ptsH5];
ptsD5 := &join[{@ D5![pt[1],pt[2]] : pt in a @} : a in $1]; ptsD5;
// (-1,-2) is the point where the map is not defined,
// (-1,-4/3) is a singularity.
//
// Check if one of these points lifts to C5.
// Note that pt @@ mC5 here will return a scheme, so we have to ask
// for its points.
[Points(pt @@ mC5) : pt in ptsD5];
// So, no point lifts.
//
// It remains to show that we really have found all points on S5 (say).
//
// So we again find the Mordell-Weil group.
J5 := Jacobian(S5);
time TwoSelmerGroupData(J5);
//
// The rank bound is 1, and by standard conjectures (which imply that
// the difference between the bound and the actual rank is even), the
// rank has to be 1.
// We can find a point of infinite order to prove this:
ptJ := J5![ptsS5[1], ptsS5[2]]; ptJ;
Order(ptJ);
//
// We need the full group J5(Q). First check the torsion.
TorsionSubgroup(J5);
// ==> no torsion.
//
// So J5(Q) is isomorphic to Z, and we need to find a generator.
// Let us look for some small points on J5.
ptsJ5 := Points(J5 : Bound := 10); ptsJ5;
//
// The following finds a basis of the group they generate.
time ReducedBasis(ptsJ5);
a,b := $1;
gen := a[1]; ht := b[1,1]; // b is a matrix!
//
// So most likely, this point generates the MW group.
// Let us prove this!
// If gen is not a generator, then it is a multiple of the generator,
// whose canonical height therefore is at most ht/4.
// Now we can get a bound on the difference between naive and canonical
// height, as follows.
bd := HeightConstant(J5 : Effort := 2); bd;
//
// This is logarithmic, so the bound for the search is
sbd := Floor(Exp(ht/4 + bd)); sbd;
//
// We need to find all points up to that naive height on J5 and check
// that they all are multiples of gen.
ptsJ5 := Points(J5 : Bound := sbd);
ptsJ5 diff {@ n*gen : n in [-10..10] @};
// OK. We have shown that J5(Q) = .
//
// Since J5(Q) is infinite, we cannot use the simple trick that worked
// for n=4.
// However, since the rank (1) is less than the genus (2), we can use
// Chabauty's method.
// There is a function available (which we will use below), but it is
// not very efficient and should be used with care. At some point, it
// should be replaced by a better, more general and more powerful version...
time Chabauty(gen, 3);
ptsS5;
// The output is to be read as follows.
// For each residue class x = infty, -2, 1 mod 3^4, there is at most one
// pair of rational points on the curve.
// Since all these possibilities are accounted for by the points we know,
// this shows that there are no other points.
//
// Therefore, we have proved another result.
// THEOREM. There are no rational 5-cycles under z |--> z^2 + c.
__