CHAPTER III ELLIPTIC CURVE ALGORITHMS

3.1 Terminology and notation

For reference in the following sections, we collect here the notation, terminology and formulae concerning elliptic curves which we will use throughout this chapter.

An elliptic curve E defined over ℚ has an equation or model of the form

(3.1.1)  E:  y^{2} + a_{1}xy + a_{3}y = x^{3} + a_{2}x^{2} + a_{4}x + a_{6}

where the coefficients a_{i} ∈ ℚ. We call such an equation a Weierstrass equation for E, and denote this model by [a_{1}, a_{2}, a_{3}, a_{4}, a_{6}]. We say that (3.1.1) is integral or defined over ℤ if all the a_{i} are in ℤ. From these coefficients we derive the auxiliary quantities

b_{2} = a_{1}^{2} + 4a_{2},
b_{4} = a_{1}a_{3} + 2a_{4},
b_{6} = a_{3}^{2} + 4a_{6},
b_{8} = a_{1}^{2}a_{6} − a_{1}a_{3}a_{4} + 4a_{2}a_{6} + a_{2}a_{3}^{2} − a_{4}^{2},

the invariants

c_{4} = b_{2}^{2} − 24b_{4},
c_{6} = −b_{2}^{3} + 36b_{2}b_{4} − 216b_{6},

the discriminant

Δ = −b_{2}^{2}b_{8} − 8b_{4}^{3} − 27b_{6}^{2} + 9b_{2}b_{4}b_{6},

and the j-invariant

j = c_{4}^{3}/Δ,

which are related by the identities

4b_{8} = b_{2}b_{6} − b_{4}^{2}  and  1728Δ = c_{4}^{3} − c_{6}^{2}.

The discriminant Δ must be non-zero for the curve defined by equation (3.1.1) to be nonsingular and hence an elliptic curve. The j-invariant is (as its name suggests) invariant under isomorphism; elliptic curves with the same j are called twists: they are isomorphic over an algebraic extension, but not necessarily over ℚ. The invariants c_{4} and c_{6} are sufficient to determine E up to isomorphism (over ℚ) since E is isomorphic to

Y^{2} = X^{3} − 27c_{4}X − 54c_{6}.

The most general isomorphism from E to a second curve E′ given by an equation of the form (3.1.1), which we usually think of as a change of coordinates on E itself, is T (r, s, t, u), given by

(3.1.2)  \begin{aligned}x & = u^{2}x′ + r \\ y & = u^{3}y′ + su^{2}x′ + t \end{aligned}

where r, s, t ∈ ℚ and u ∈ ℚ^{∗}. The effect of T (r, s, t, u) on the coefficients a_{i} is given by

(3.1.3)  \begin{aligned}ua′_{1} & = a_{1} + 2s \\ u^{2}a′_{2} & = a_{2} − sa_{1} + 3r − s^{2} \\ u^{3}a′_{3} & = a_{3} + ra_{1} + 2t \\ u^{4}a′_{4} & = a_{4} − sa_{3} + 2ra_{2} − (t + rs)a_{1} + 3r^{2} − 2st \\ u^{6}a′_{6} & = a_{6} + ra_{4} + r^{2}a_{2} + r^{3} − ta_{3} − t^{2} − rta_{1} \end{aligned}

so that

u^{4}c′_{4} = c_{4},  u^{6}c′_{6} = c_{6},  u^{12}Δ′ = Δ and j′ = j.

The transformations T (0, 0, 0, u) we will refer to as scaling transformations; these have the effect of dividing each coefficient a_{i} by u^{i}, and similarly for each of the other quantities, according to its weight. Here a_{i}, b_{i} and c_{i} have weight i, while Δ has weight 12 and j has weight 0. By applying T (0, 0, 0, u) for suitable u we can always transform to an integral model; all the invariants are then integral, except (possibly) for j. Among such integral models, those for which the positive integer |Δ| is minimal are called global minimal models for E. We will give in the next section a simple algorithm for finding such a model, given the invariants c_{4} and c_{6} of any model. Clearly, isomorphisms between minimal models must have u = ±1 and r, s, t ∈ ℤ. We may normalize so that a_{1}, a_{3} ∈ {0, 1} and a_{2} ∈ {−1, 0, 1}, by suitable choice of s, r and t (in that order), as may be seen from (3.1.3). Such an equation will be called reduced, and it is not hard to show that it is unique: the only transformation other than the identity T (0, 0, 0, 1) from a reduced model to any another reduced model is the transformation T (0, −a_{1}, −a_{3}, −1), which takes any model to itself; this is just the negation map (x, y) ↦ (x, y − a_{1}x − a_{3}) from the curve to itself. Thus every elliptic curve E defined over ℚ has a unique reduced minimal model. This fact makes it very easy to recognize curves: in Table 1 we give the coefficients of such a model for each of the curves there.

Given integers c_{4} and c_{6}, two questions arise: is there a curve over ℚ with these invariants, and is it minimal? Clearly we must have c_{4}^{3} − c_{6}^{2} = 1728Δ with Δ ≠ 0. A solution to the first problem is given by Kraus in Proposition 2 of [30], which states the following.

Proposition 3.1.1. Let c_{4}, c_{6} be integers such that Δ = (c_{4}^{3} − c_{6}^{2})/1728 is a non-zero integer. In order for there to exist an elliptic curve E with a model (3.1.1) defined over ℤ having invariants c_{4} and c_{6}, it is necessary and sufficient that

(1) c_{6} ≢ ±9 (mod 27);

(2) either c_{6} ≡ −1 (mod 4), or c_{4} ≡ 0 (mod 16) and c_{6} ≡ 0, 8 (mod 32).

The conditions of Proposition 3.1.1 will be referred to as Kraus’s conditions. If we are given integers c_{4} and c_{6} satisfying these conditions, we can recover the coefficients a_{i} of the reduced model of the curve with c_{4} and c_{6} as invariants, using the formulae already given in Chapter 2, Section 14, which we repeat here for convenience:

b_{2} = −c_{6} mod 12 ∈ {−5, …, 6};
b_{4} = (b_{2}^{2} − c_{4})/24;
b_{6} = (−b_{2}^{3} + 36b_{2}b_{4} − c_{6})/216;
a_{1} = b_{2} mod 2 ∈ {0, 1};
a_{3} = b_{6} mod 2 ∈ {0, 1};
a_{2} = (b_{2} − a_{1})/4;
a_{4} = (b_{4} − a_{1}a_{3})/2;
a_{6} = (b_{6} − a_{3})/4.

To see this, we may assume that we are seeking coefficients of a reduced model; then b_{2} ∈ {−4, −3, 0, 1, 4, 5}, and we have −c_{6} ≡ b_{2}^{3} ≡ b_{2} (mod 12). The rest is easy; provided that c_{4} and c_{6} satisfy Kraus’s conditions, all the divisions will be exact.

In the following section we answer the second question by giving an algorithm for computing the reduced coefficients of a minimal model for any curve E, given either integral invariants satisfying Kraus’s conditions, or any integral model for E. We simply determine the maximal integer u such that c′_{4} = c_{4}/u^{4} and c′_{6} = c_{6}/u^{6} satisfy Kraus’s conditions, and then compute the reduced coefficients a′_{i} from these. As with many questions concerning elliptic curves, most of the work goes into determining the powers of 2 and 3 which divide u.

We will assume without further discussion that on any given curve E, points may be added and multiples taken, using standard formulae. The Mordell–Weil group of all rational points on E will be denoted E(ℚ) as usual. If n is a positive integer, we denote by E(ℚ)[n] the subgroup of rational points of order dividing n, which is the kernel of the multiplication map from E to itself.

3.2 The Kraus–Laska–Connell algorithm and Tate’s algorithm

In this section we give two algorithms. The first was originally given by Laska in [34], and finds a minimal model for a curve E, starting from an integral equation. Essentially the algorithm was to test all positive integers u such that u^{−4}c_{4} and u^{−6}c_{6} are integral, to see if they are the invariants of a curve defined over ℤ. Using Kraus’s conditions (see Proposition 3.1.1 above), this procedure can be simplified, since it is possible to compute in advance the exponent d_{p} of each prime p in the minimal discriminant, and hence compute u at the start. The usual formulae then give the coefficients a_{i} of the reduced model. Our formulation of the resulting algorithm over ℤ is similar to that given in [10], where more general rings are considered: in particular an explicit algorithm is given there for finding local minimal models over arbitrary number fields, and hence global minimal models where they exist. Over ℤ, the algorithm is extremely simple.

In the pseudocode below,

ord(p,n) gives the power of the prime p which divides the non-zero integer n; floor(x) gives the integral part of the real number x; a mod p gives the residue of a modulo p lying in the range −\frac{1}{2}p < a ≤ \frac{1}{2}p; in particular,

when p = 2 or 3 this gives a residue in {0, 1} or {−1, 0, 1} respectively. Also inv(a,p) gives the inverse of a modulo p, assuming that gcd(a,p)=1.

The Laska–Kraus–Connell Algorithm

INPUT:        c4, c6 (integer invariants of an elliptic curve E).
OUTPUT:       a1, a2, a3, a4, a6 (coefficients of a reduced minimal model for E).
 1.  BEGIN
 2. Δ=(c43-c62)/1728;

(Compute scaling factor u)

 3.  u = 1; g = gcd(c62,Δ);
 4.  p list = prime divisors(g);
 5.  FOR p IN p list DO
 6.  BEGIN
 7.       d = floor(ord(p,g)/12);
 8.       IF p=2 THEN
 9.           a = c4/2(4*d) mod 16; b = c6/2(6*d) mod 32;
10.           IF (b mod 4 ≠ -1) AND NOT (a=0 AND (b=0 OR b=8))

11.           THEN d = d-1
12.           FI
13.      ELIF p=3 THEN IF ord(3,c6)=6*d+2 THEN d = d-1 FI
14.      FI;
15.      u = u*pd
16.  END;

(Compute minimal equation)

17.  c4 = c4/u4; c6 = c6/u6;
18.  b2 = -c6 mod 12; b4 = (b22-c4)/24; b6 = (-b23+36*b2*b4-c6)/216;
19.  a1 = b2 mod 2;
20.  a3 = b6 mod 2;
21.  a2 = (b2-a1)/4;
22.  a4 = (b4-a1*a3)/2;
23.  a6 = (b6-a3)/4
24.  END

Next we turn to Tate’s algorithm itself. The standard reference for this is Tate’s ‘letter to Cassels’ [65], which appeared in the Antwerp IV volume [2]. There is also a full account in the second volume of Silverman’s book [61, Section IV.9]. It may be applied to an integral model of a curve E and a prime p, to give the following data:

• The exponent f_{p} of p in the conductor N of E (see below);

• the Kodaira symbol of E at p, which classifies the type of reduction of E at p (see [47] or [61, Section IV.9]); these are: I_{0} for good reduction; I_{n} (n > 0) for bad multiplicative reduction; and types I_{n}^{∗}, II, III, IV, II^{∗}, III^{∗} and IV^{∗} for bad additive reduction.

• the local index c_{p} = [E(ℚ_{p}) : E^{0}(ℚ_{p})], where E^{0}(ℚ_{p}) is the subgroup of the group E(ℚ_{p}) of p-adic points of E, consisting of those points whose reduction modulo p is non-singular. (That this index is finite is implied by the correctness of the algorithm, as observed by Tate in [65].)

In addition, the algorithm detects whether the given model is non-minimal at p, and if so, returns a model which is minimal at p. Thus by applying it in succession with all the primes dividing the discriminant of the original model, one can compute a minimal model at the same time as computing the conductor and the other local reduction data. In practice this makes the Laska–Kraus–Connell algorithm redundant, though much simpler to implement and use if all one needs is the standard model for a curve E.

The conductor N of an elliptic curve E defined over ℚ is defined to be

N = ∏_{p} p^{f_{p}}

where f_{p} = ord_{p}(Δ) + 1 − n_{p} and n_{p} is the number of irreducible components on the special fibre of the minimal Néron model of E at p. This Néron model is a more sophisticated object than we wish to discuss here (see [47] or [61] for details): one has to consider E as a scheme over Spec(ℤ_{p}), and then resolve the singularity at p, to obtain a scheme whose generic fibre is E/ℚ_{p} and whose special fibre is a union of curves over ℤ/pℤ. In terms of a minimal model for E over ℤ, all may be computed very simply except when p = 2 or p = 3 as follows:

f_{p} = 0 if p ∤ Δ;

f_{p} = 1 if p | Δ and p ∤ c_{4} (then n_{p} = ord_{p}(Δ));

f_{p} ≥ 2 if p | Δ and p | c_{4}; moreover, f_{p} = 2 in this case when p ≠ 2, 3.

To obtain the value of f_{p} in the remaining cases, and to obtain the Kodaira symbol and the local index c_{p}, we use Tate’s algorithm itself.

In [65], the algorithm is given for curves defined over an arbitrary discrete valuation ring. To apply it to a curve defined over the ring of integers R of a number field K at a prime ideal 𝔭, one would in general have to work in the localization of R at 𝔭; here we can work entirely over ℤ, since ℤ is a principal ideal domain. We have added to the presentation in [65] the explicit coordinate transformations T (r, s, t, u) which are required during the course of the algorithm to achieve divisibility of the coefficients a_{i} by various power of p. In practice one would ignore the transformations which had taken place while processing each p, unless a scaling by p had taken place on discovering that the model was non-minimal. The most complicated part of the algorithm is the branch for reduction type I_{m}^{∗}, where one successively refines the model p-adically until certain auxiliary quadratics have distinct roots modulo p. This requires careful book-keeping. The presentation given here closely follows our own implementation of the algorithm, which in turn owes much to an earlier Fortran program written by Pinch. The following sub-procedures are used:

compute invariants computes the b_{i}, c_{i} and Δ from the coefficients a_{i}. Note that c_{4}, c_{6} and Δ do not change unless a scaling is required, since all other transformations have u = 1.

transcoord(r,s,t,u) applies the coordinate transformation formulae of the previous section to obtain new values for the a_{i} and other quantities. All calls to this procedure have u = 1 except when rescaling a non-minimal equation. In each case we first compute suitable values of r, s and t; usually this requires a separate branch if p = 2 or p = 3.

quadroots(a,b,c,p) returns TRUE if the quadratic congruence ax^{2} + bx + c ≡ 0 (mod p) has a solution, and FALSE otherwise. This is used in determining the value of the index c_{p}.

nrootscubic(b,c,d,p) returns the number of roots of the cubic congruence x^{3} +bx^{2} +cx+ d ≡ 0 (mod p).

Tate’s Algorithm

INPUT:        a1, a2, a3, a4, a6 (integer coefficients of E); p (prime).
OUTPUT:       Kp   (Kodaira symbol)
              fp   (Exponent of p in conductor)
              cp   (Local index)
 1.  BEGIN
 2.  compute invariants(b2,b4,b6,b8,c4,c6,Δ);
 3.  n = ord(p,Δ);

(Test for type I_{0})

 4.  IF n=0 THEN Kp = "I0"; fp = 0; cp = 1; EXIT FI;

(Change coordinates so that p | a_{3}, a_{4}, a_{6})

 5.  IF p=2 THEN
 6.       IF p|b2
 7.       THEN r = a4 mod p; t = r*(1+a2+a4)+a6 mod p
 8.       ELSE r = a3 mod p; t = r+a4 mod p
 9.       FI
10.  ELIF p=3 THEN
11.      IF p|b2 THEN r = -b6 mod p ELSE r = -b2*b4 mod p FI;
12.      t = a1*r+a3 mod p
13.  ELSE
14.      IF p|c4 THEN r = -inv(12,p)*b2 ELSE r = -inv(12*c4,p)*(c6+b2*c4) FI;
15.      t = -inv(2,p)*(a1*r+a3);
16.      r = r mod p; t = t mod p
17.  FI;

18.  transcoord(r,0,t,1);

(Test for types I_{n}, II, III, IV)

19.  IF p∤c4 THEN
20.      IF quadroots(1,a1,-a2,p) THEN cp = n ELIF 2|n THEN cp = 2 ELSE cp = 1 FI;
21.      Kp = "In"; fp = 1; EXIT
22.  FI;
23.  IF p2∤a6 THEN Kp = "II"; fp = n; cp = 1; EXIT;
24.  IF p3∤b8 THEN Kp = "III"; fp = n-1; cp = 2; EXIT;
25.  IF p3∤b6 THEN
26.      IF quadroots(1,a3/p,-a6/p2,p) THEN cp = 3 ELSE cp = 1 FI;
27.      Kp = "IV"; fp = n-2; EXIT
28.  FI;

(Change coordinates so that p | a_{1}, a_{2}; p^{2} | a_{3}, a_{4}; p^{3} | a_{6})

29.  IF p=2
30.  THEN s = a2 mod 2; t = 2*(a6/4 mod 2)
31.  ELSE s = -a1*inv(2,p); t = -a3*inv(2,p)
32.  FI;
33.  transcoord(0,s,t,1);

(Set up auxiliary cubic T^{3} + bT^{2} + cT + d)

34.  b = a2/p; c = a4/p2; d = a6/p3;
35.  w = 27*d2-b2*c2+4*b3*d-18*b*c*d+4*c3;
36.  x = 3*c-b2;

(Test for distinct roots: type I_{0}^{∗})

37.  IF p∤w THEN Kp = "I*0"; fp = n-4; cp = 1+nrootscubic(b,c,d,p); EXIT

(Test for double root: type I_{m}^{∗})

38.  ELIF p∤x THEN

(Change coordinates so that the double root is T ≡ 0)

39.      IF p=2 THEN r = c ELIF p=3 THEN r = b*c ELSE r = (b*c-9*d)*inv(2*x,p) FI;
40.      r = p*(r mod p);
41.      transcoord(r,0,0,1);

(Make a_{3}, a_{4}, a_{6} repeatedly more divisible by p)

42.      m = 1; mx = p2; my = p2; cp = 0;
43.      WHILE cp=0 DO
44.      BEGIN
45.           xa2 = a2/p; xa3 = a3/my; xa4 = a4/(p*mx); xa6 = a6/(mx*my);
46.           IF p∤(xa32+4*xa6) THEN
47.               IF quadroots(1,xa3,-xa6,p) THEN cp = 4 ELSE cp = 2 FI
48.           ELSE
49.               IF p=2 THEN t = my*xa6 ELSE t = my*((-xa3*inv(2,p)) mod p) FI;
50.               transcoord(0,0,t,1);
51.               my = my*p; m = m+1;
52.               xa2 = a2/p; xa3 = a3/my; xa4 = a4/(p*mx); xa6 = a6/(mx*my);
53.               IF p∤(xa42-4*xa2*xa6) THEN
54.                    IF quadroots(xa2,xa4,xa6,p) THEN cp = 4 ELSE cp = 2 FI
55.               ELSE
56.                    IF p=2 THEN r = mx*(xa6*xa2 mod 2)

57.                    ELSE r = mx*(-xa4*inv(2*xa2,p) mod p)
58.                    FI;
59.                    transcoord(r,0,0,1);
60.                    mx = mx*p; m = m+1
61.               FI
62.           FI
63.      END;
64.      fp = n-m-4; Kp = "I*m"; EXIT
65.  ELSE

(Triple root case: types II^{∗}, III^{∗}, IV^{∗} or non-minimal) (Change coordinates so that the triple root is T ≡ 0)

66.      IF p=3 THEN rp = -d ELSE rp = -b*inv(3,p) FI;
67.      r = p*(rp mod p);
68.      transcoord(r,0,0,1);
69.      x3 = a3/p2; x6 = a6/p4;

(Test for type IV^{∗})

70.      IF p∤(x32+4*x6) THEN
71.           IF quadroots(1,x3,-x6,p) THEN cp = 3 ELSE cp = 1 FI;
72.           Kp = "IV*"; fp = n-6; EXIT
73.      ELSE

(Change coordinates so that p^{3} | a_{3}, p^{5} | a_{6})

74.           IF p=2 THEN t = x6 ELSE t = x3*inv(2,p) FI;
75.           t = -p2*(t mod p);
76.           transcoord(0,0,t,1);

(Test for types III^{∗}, II^{∗})

77.           IF p4∤a4 THEN Kp = "III*"; fp = n-7; cp = 2; EXIT
78.           ELIF p6∤a6 THEN Kp = "II*"; fp = n-8; cp = 1; EXIT
79.           ELSE

(Equation non-minimal: divide each a_{i} by p^{i} and start again)

80.               transcoord(0,0,0,p); restart
81.           FI
82.  FI
83.  END

In Table 1 we will give the local reduction data for each curve at each ‘bad’ prime (dividing the discriminant of the minimal model). We also give the factorization of the minimal discriminant and of the denominator of j, as in the earlier tables. To save space we omit the c_{4} and c_{6} invariants, which are easily computable from the coefficients a_{i}.

3.3 Computing the Mordell–Weil group I: finding torsion points

In this and the next three sections we will discuss the question of determining the Mordell– Weil group E(ℚ) of rational points on an elliptic curve E defined over ℚ. This group is finitely generated, by Mordell’s Theorem, and hence has the structure

E(ℚ) = T × F

where T is the finite torsion subgroup E(ℚ)_{tors} of E(ℚ) consisting of the points of finite order, and F is free abelian of some rank r ≥ 0:

F ≅ ℤ^{r}.

The problem of computing E(ℚ) thus subdivides into several parts:

• computing the torsion T;

• computing the rank r;

• finding r independent points of infinite order;

• computing a ℤ-basis for the free part F.

A related task is to compute the regulator R(E(ℚ)) (defined below); for this and for the latter two steps we will also need to compute the canonical height ĥ(P) of points P ∈ E(ℚ), and hence the height pairing ĥ(P, Q).

In this section we will treat the easiest of these problems, that of finding the torsion points. In fact, these can be found as a byproduct of the more general search for points on the curve, since their naive height can be bounded (see the remark before Lemma 3.5.2). However, it is also useful to have a self-contained method for determining the torsion.

Using the fact that E(ℝ) is isomorphic either to the circle group S^{1} (when Δ < 0) or to S^{1} × C_{2} (when Δ > 0), where C_{k} denotes a cyclic group of order k, together with the fact that all finite subgroups of S^{1} are cyclic, we see that T is isomorphic either to C_{k} or to C_{2k} × C_{2} for some k ≥ 1, the latter only being possible when Δ is positive. The number of possible values of k is finite: by a theorem of Mazur [39],[40], a complete list of possible structures of T is

C_{k}  for  1 ≤ k ≤ 10 or k = 12;
C_{2k} × C_{2}  for  1 ≤ k ≤ 4.

To determine the torsion subgroup of an elliptic curve defined over ℚ, we may use a form of the Lutz–Nagell Theorem. (The situation is more complicated over number fields other than ℚ, on account of the ramified primes.) The first step is to find a model for the curve in which all torsion points are integral. For this it suffices to complete the square (if necessary) to eliminate the xy and y terms, at the expense of a scaling by u = 2. Then for P = (x, y) a torsion point, we can use the fact that both P and 2P are integral to bound y. For the first step, the following result may be found in [33, Section III.1] and [28, Theorem 5.1]. The original form of this result, due independently to Lutz [36] and Nagell [46], was for curves of the form y^{2} = x^{3} + ax + b, with no x^{2} term. While such an equation may be obtained by completing the cube, this would involve a further scaling of coordinates, and so would lead to larger numbers. If a_{1} = a_{3} = 0 we can apply the following result directly; otherwise, put a = b_{2}, b = 8b_{4} and c = 16b_{6}.

Proposition 3.3.1. Let E be an elliptic curve defined over ℚ, given by an equation

(3.3.1)  y^{2} = f(x) = x^{3} + ax^{2} + bx + c

where a, b, c ∈ ℤ. If P = (x, y) ∈ E(ℚ) has finite order, then x, y ∈ ℤ.

Next we bound the y coordinate of a torsion point P = (x, y) (see [33, Theorem 1.4]).

Proposition 3.3.2. Let E be as in (3.3.1). If P = (x_{1}, y_{1}) has finite order in E(ℚ) then either y_{1} = 0 or y_{1}^{2} | Δ_{0}, where

Δ_{0} = 27c^{2} + 4a^{3}c + 4b^{3} − a^{2}b^{2} − 18abc.

Proof. If 2P = 0 then y_{1} = 0, since −P = (x_{1}, −y_{1}). Otherwise 2P = (x_{2}, y_{2}) with x_{2}, y_{2} ∈ ℤ by Proposition 3.3.1. Using the addition formula on E we find that 2x_{1} + x_{2} = m^{2} − a where m = f′(x_{1})/2y_{1} is the slope of the tangent to E at P. Hence m ∈ ℤ, so that y_{1} | f′(x_{1}). Using y_{1}^{2} = f(x_{1}), this implies that y_{1}^{2}|Δ_{0}, since

Δ_{0} = (−27f(x) + 54c + 4a^{3} − 18ab)f(x) + (f′(x) + 3b − a^{2})f′(x)^{2}. □

This gives us a finite number of values of y to check; for each, we attempt to solve the cubic for x ∈ ℤ, to obtain all torsion points on E. Note that we are actually determining all points P such that both P and 2P are integral (in the possibly scaled model for E), which includes all torsion points, but may also include points of infinite order. To determine whether a given integral point has finite or infinite order, we simply compute multiples mP successively until either mP = 0, in which case P has order m, or mP is not integral, in which case P has infinite order. This does not take long, as the maximum possible order for a torsion point is 12 by Mazur’s theorem. If we find points of infinite order at this stage we keep a note of them for later use (see Section 3.5).

The quantity Δ_{0} is related to the discriminant Δ of the curve (3.3.1) by Δ = −16Δ_{0}. If this is large, there may be many values of y_{0} to check when we apply the preceding Proposition to determine the torsion on a given curve. It is possible to save time by using a further result, which states that for an odd prime p of good reduction (that is, p ∤ 2Δ), the reduction map from E(ℚ)_{tors} to E(ℤ/pℤ) is injective. For more details, and worked examples, see either [58, Section VIII.7] or [28, Section V.1].

If we want to know the structure of T and not just its order, note that from Mazur’s theorem the only ambiguous cases are when T has order 4k = 4, 8 or 12 and Δ > 0; we can always tell apart the groups C_{4k} and C_{2} × C_{2k} as the former has only one element of order 2 while the latter has three, and this number is the number of rational (integer) roots of f(x).

To solve the cubic equations f(x) = y^{2} for x, given y, we use the classical formula of Cardano (see any algebra textbook) to find the complex roots (which we also need in computing the periods in section 3.7 below), and if any of these are real and close to integers we check them using exact integer arithmetic. Testing all divisors of the constant term can be too timeconsuming, as it involves factorization of the numbers y^{2} − c which may be very large.

Here is the algorithm in pseudocode; for simplicity we only give it for curves with no xy or y terms; in the general case, one works internally with points on a scaled model (including the calculation of the order), converting back to the original model on output. Since we know in advance that no point will have order greater than 12, when computing the order of a point we simply use repeated addition until we reach a non-integral point or the identity O. The subroutine order(P) returns 0 for a point of infinite order. Also: square part(Δ) returns the largest integer whose square divides Δ; integer roots returns a list of the integer roots of a cubic with integral coefficients; and integral(x) tests whether its (rational) argument is integral.

Algorithm for finding all torsion points

INPUT:        a,b,c (integer coefficients of a nonsingular cubic).
OUTPUT:       A list of all torsion points on y2=x3+ax2+bx+c, with orders.
 1.  BEGIN
 2. Δ=27*c2+4*a3*c+4*b3-a2*b2-18*a*b*c;
 3.  y list=positive divisors(square part(Δ)) ∪ {0};
 4.  FOR y IN y list DO
 5.  BEGIN

 6.       x list=integer roots(x3+a*x2+b*x+c-y2);
 7.       FOR x IN x list DO
 8.       BEGIN
 9.           P=point(x,y);
10.           n=order(P);
11.           IF n>0 THEN OUTPUT P,n FI
12.      END
13.  END
14.  END

(Subroutine to compute order of a point)

SUBROUTINE order(P)

 1.  BEGIN
 2.  n=1; Q=P;
 3.  WHILE integral(x(Q)) AND Q≠O DO
 4.  BEGIN
 5.       n = n+1; Q = Q+P
 6.  END;
 7.  IF Q≠O THEN n=0 FI;
 8.  RETURN n
 9.  END

3.4 Heights and the height pairing

In this section we will show how to compute the canonical height ĥ(P) of a point P ∈ E(ℚ), and hence the height pairing

ĥ(P, Q) = \frac{1}{2}(ĥ(P + Q) − ĥ(P) − ĥ(Q)).

We will use this in the following section to find dependence relations among finite sets of points of infinite order, when we are computing a ℤ-basis {P_{1}, …, P_{r}} for the free abelian group E(ℚ)/T. Also, the regulator R(E) is given by the determinant of the height pairing matrix:

R(E) = |det(ĥ(P_{i}, P_{j}))| .

The canonical height ĥ is a real-valued quadratic form on E(ℚ). It differs by a bounded amount (with a bound dependent on E but not on the point P) from the naive or Weil height h(P). For a point P = (x, y) = (a/c^{2}, b/c^{3}) ∈ E(ℚ) with a, b, c ∈ ℤ and gcd(a, c) = 1 = gcd(b, c), the latter is defined to be

h(P) = log max{|a|, c^{2}}.

Now the canonical height may be defined as ĥ(P) = lim_{n→∞} 4^{−n}h(2^{n}P), but this is not practical for computational purposes. For the theory of heights on elliptic curves, see [58, Chapter VIII]. Later (in the next section) we will need an explicit bound on the difference between ĥ(P) and h(P).

The height algorithms in this section are taken from Silverman’s paper [59]. The global height ĥ(P) is defined as a sum of local heights:

(3.4.1)  ĥ(P) = ∑_{p≤∞} ĥ_{p}(P).

Here the sum is over all finite primes p and the ‘infinite prime’ ∞ coming from the real embedding of ℚ. (Over a general number field, there would in general be several of these infinite primes, including complex ones, and the local heights need to be multiplied by certain multiplicities: see [59]).

A remark about normalization^{1}: the canonical height must be suitably normalized. In the literature there are two normalizations used, one of which is double the other and is the one appropriate for the Birch–Swinnerton-Dyer conjecture (resulting in a regulator 2^{r} times as large). In Silverman’s paper he uses the other (smaller) normalization. Thus all the formulae here are double those in the paper [59].

The following proposition, which is Theorem 5.2 of [59] (for curves over general number fields) specialized to the case of a curve defined over ℚ, also applies to a curve defined over ℚ_{p} and to a point P = (x, y) ∈ E(ℚ_{p}). In the proposition, we refer to the functions ψ_{2} and ψ_{3} defined on E by

ψ_{2}(P) = 2y + a_{1}x + a_{3},  and  ψ_{3}(P) = 3x^{4} + b_{2}x^{3} + 3b_{4}x^{2} + 3b_{6}x + b_{8};

thus, ψ_{2} vanishes at the 2-torsion points of E and ψ_{3} at the 3-torsion.

Proposition 3.4.1. Let E be an elliptic curve defined over ℚ given by a standard Weierstrass equation (3.1.1) which is minimal at p, and let P = (x, y) ∈ E(ℚ).

(a) If

ord_{p}(3x^{2} + 2a_{2}x + a_{4} − a_{1}y) ≤ 0  or  ord_{p}(2y + a_{1}x + a_{3}) ≤ 0

then

ĥ_{p}(P) = max{0, −ord_{p}(x)} log p.

(b) Otherwise, if ord_{p}(c_{4}) = 0 then set N = ord_{p}(Δ) and M = min\{ord_{p}(ψ_{2}(P), \frac{1}{2}N\}; then

ĥ_{p}(P) = \frac{M(M − N)}{N} log p.

(c) Otherwise, if ord_{p}(ψ_{3}(P)) ≥ 3ord_{p}(ψ_{2}(P)) then

ĥ_{p}(P) = −\frac{2}{3}ord_{p}(ψ_{2}(P)) log p.

(d) Otherwise

ĥ_{p}(P) = −\frac{1}{4}ord_{p}(ψ_{3}(P)) log p.

The first case in Proposition 3.4.1 covers primes p where the point P has good reduction (including all primes where E has good reduction, as well as those where the reduced curve is singular but P does not reduce to the singular point). In the other three cases, P has singular reduction, and the reduction of E at p is multiplicative, additive of types IV or IV^{∗}, and additive of types III, III^{∗} and I_{m}^{∗} respectively.

Hence for each point P, the local height ĥ_{p}(P) = 0 if p divides neither the discriminant Δ nor c, where c^{2} is the denominator of the x-coordinate of the point P. In all cases, ĥ_{p}(P) is a rational multiple of log(p). The total contribution from the primes dividing c in the global height ĥ(P) is therefore (from case (a) of the Proposition) simply 2 log(c), and we have the following formula, better for practical computation than (3.4.1) since we do not have to factorize c:

1I am grateful to Gross for explaining this to me, after I found that apparently the two sides of the Birch– Swinnerton-Dyer conjecture disagreed by a factor of 2^{r}!

(3.4.2)  ĥ(P) = ĥ_{∞}(P) + 2 log(c) + ∑_{p|Δ,p∤c} ĥ_{p}(P).

This formula appears in [62], where it is shown how to compute ĥ(P) using little (or no) factorization of Δ, which can be useful in certain situations. We refer the reader to [62] for details.

An algorithm for computing the local height at a finite prime p is given by the following:

Silverman’s algorithm for computing local heights: finite primes

INPUT:   a1, a2, a3, a4, a6 (integer coefficients of a minimal model for E).

         x,y (rational coordinates of a point P on E).
         p (a prime).

OUTPUT:  the local height of P at p.
 1.  BEGIN
 2.  compute invariants(b2,b4,b6,b8,c4,Δ);
 3.  N = ord(p,Δ);
 4.  A = ord(p,3*x2+2*a2*x+a4-a1*y);
 5.  B = ord(p,2*y+a1*x+a3);
 6.  C = ord(p,3*x4+b2*x3+3*b4*x2+3*b6*x+b8);
 7.  M = min(B,N/2);
 8.  IF A ≤ 0 OR B ≤ 0 THEN L = max(0,-ord(p,x))
 9.  ELSE IF ord(p,c4)=0 THEN L = M*(M-N)/N
10.  ELSE IF C ≥ 3*B THEN L = -2*B/3
11.  ELSE L = -C/4
12.  FI;
13.  RETURN L*log(p)
14.  END

We must also compute the local component of the height at the infinite prime, ĥ_{∞}(P). The method here originated with Tate, but was amended by Silverman in [59] to improve convergence, and to apply also to complex valuations. Tate in [66] expressed ĥ_{∞}(P) as a series

ĥ_{∞}(P) = log |x| + \frac{1}{4} ∑_{n=0}^{∞} 4^{−n}c_{n}

where the coefficients c_{n} are bounded provided that no point on E(ℝ) has x-coordinate zero. Of course, over ℝ one can shift coordinates to ensure that this condition holds, but the resulting series can have poor convergence properties, and this trick will not work over ℂ. Silverman’s solution is to use alternately the parameters x and x′ = x + 1, switching between them (and between the two associated series c_{n} and c′_{n}) whenever |x| or |x′| becomes small (less than 1/2). The series of coefficients c_{n} is obtained by repeated doubling of the point P, working with t = 1/x or t′ = 1/x′ as local parameter. The result is a new series of the above type in which the error in truncating before the Nth term is O(4^{−N}), with an explicit constant. In fact (see [59, Theorem 4.2]) the error is less than \frac{1}{2}10^{−d}, giving a result correct to d decimal places, if

N ≥ \frac{5}{3}d + \frac{1}{2} + \frac{3}{4} log(7 + \frac{4}{3} log H + \frac{1}{3} log max\{1, |Δ|^{−1}\})

where

H = max{4, |b_{2}|, 2|b_{4}|, 2|b_{6}|, |b_{8}|}.

The last term vanishes for curves defined over ℤ, since then we have |Δ| > 1.

In the algorithm which we now give, the quantities b2', b4', b6' and b8' are those associated with the shifted model of E with x′ = x+1; the switching flag beta indicates which model we are currently working on; mu holds the current partial sum; f holds the negative power of 4.

Silverman’s algorithm for computing local heights: real component

INPUT:   a1, a2, a3, a4, a6 (integer coefficients of a minimal model for E).

         x (x-coordinate of a point P on E).
         d (number of decimal places required).

OUTPUT:  the real local height of P.
 1.  BEGIN
 2.  compute invariants(b2,b4,b6,b8);
 3.  H = max(4,|b2|,2*|b4|,2*|b6|,|b8|);
 4.  b2' = b2-12; b4' = b4-b2+6; b6' = b6-2*b4+b2-4; b8' = b8-3*b6+3*b4-b2+3;
 5.  N = ceiling((5/3)*d + (1/2) + (3/4)*log(7+(4/3)*log(H)));
 6.  IF |x|<0.5 THEN t = 1/(x+1); beta = 0 ELSE t = 1/x; beta = 1 FI;
 7.  mu = -log|t|; f = 1;
 8.  FOR n = 0 TO N DO
 9.  BEGIN
10.      f = f/4;
11.      IF beta=1 THEN
12.           w = b6*t4+2*b4*t3+b2*t2+4*t;
13.           z = 1-b4*t2-2*b6*t3-b8*t4;
14.           zw = z+w
15.      ELSE
16.           w = b6'*t4+2*b4'*t3+b2'*t2+4*t;
17.           z = 1-b4'*t2-2*b6'*t3-b8'*t4;
18.           zw = z-w
19.      FI;
20.      IF |w| ≤ 2*|z|
21.      THEN mu = mu+f*log|z|; t = w/z
22.      ELSE mu = mu+f*log|zw|; t = w/zw; beta = 1-beta
23.      FI
24.  END;
25.  RETURN mu
26.  END

Finally, to compute the global height ĥ(P), we simply add to the infinite local height ĥ_{∞}(P) the finite local heights ĥ_{p}(P) for all primes p dividing either Δ or the denominator of x(P). Using (3.4.2) this leads to the following algorithm.

Algorithm for computing global canonical heights

INPUT:   a1, a2, a3, a4, a6 (integer coefficients of a minimal model for E).

         P=(x,y) (a rational point P on E).

OUTPUT:  the global canonical height ĥ(P) of P.
 1.  BEGIN
 2. Δ= discr(a1,a2,a3,a4,a6);
 3.  d = denom(x);
 4.  h = real height(P) + log(d);
 5.  p list = prime divisors(Δ);
 6.  FOR p IN p list DO
 7.  BEGIN
 8.       IF p∤d THEN h = h + local height(p,P) FI
 9.  END;
10.  RETURN h
11.  END

3.5 The Mordell–Weil group II: generators

In this section we will show how we look for rational points of infinite order on an elliptic curve E. In compiling the tables, we usually knew the rank r in advance so that we knew how many independent points to expect to find (and only looked for such points when we knew that r > 0); however, this procedure is also useful as an open-ended search when we do not know the rank, as obviously it can provide us with a lower bound for r.

The procedure divides into two parts. First, we have a searching routine which looks for points up to some bound on the naive height (equivalently, some bound on the numerator and denominator of the x-coordinate). As this routine finds points, it gives them to the second routine, which has at each stage a ℤ-basis for a subgroup A of E(ℚ)/T: initially A = 0. This second routine uses the height pairing to determine one of three possibilities: the new point P may be independent of those already found and can then be added to our cumulative list of independent points; the rank of A is thus increased by 1. Secondly, P may be an integral combination of the current basis (modulo torsion) and can then be ignored. Finally, if a multiple kP of P is an integral combination of the current basis for some k > 1, we can find a basis for a new subgroup A which contains the old A with index k. Even when we know the rank r in advance, we do not stop as soon as we have a subgroup A of rank r, since A might still have finite index in E(ℚ)/T. To close this final gap we use explicit bounds for the difference between the naive and canonical heights, such as Silverman’s result (Proposition 3.5.1) below.

The algorithm we use for the second procedure is a very general one, which can be used in many other similar situations; for example, as part of an algorithm for finding the unit group of a number field, where the first routine somehow finds units. Our algorithm is essentially the same as the ‘Algorithm for enlarging sublattices’ in the book by Pohst and Zassenhaus [50, Chapter 3.3].

A rational point P on E (given by a standard Weierstrass equation) may be written uniquely as P = (x, y) = (a/c^{2}, b/c^{3}) with integers a, b, and c satisfying gcd(a, c) = gcd(b, c) = 1 and c ≥ 1. The naive or Weil height of P is h(P) = log max{|a|, c^{2}}. Initially, we find the point of order 2 in E(ℝ) with minimal x-coordinate x_{0}; this gives a lower bound for the x-coordinates of all real points on E. We then search for points P with naive height up to some bound B by looping through positive integers c ≤ exp(B/2) and through a coprime to c in the range^{2} max{c^{2}x_{0}, − exp(B)} ≤ a ≤ exp(B). Given a and c, we attempt to solve the appropriate quadratic equation for b ∈ ℤ. To speed up this procedure, we use a quadratic sieve: for each denominator c we precompute for about 10 auxiliary sieving primes p the residue classes modulo p to which a must belong if the equation for b is to be soluble modulo p. Each candidate value of a can then first be checked to see if it is admissible modulo each sieving prime before the more time-consuming step of attempting to solve for b. This improvement to the search results in a major time saving in most cases, though for most of the curves in our tables on which we expected to find points of infinite order, such a point was found very quickly anyway. (In some cases we had already found such a point during the search for torsion points.) In practice it may be better to use composite moduli for the sieving.

2If E(ℝ) has three points of order two, with x-coordinates x_{0} < x_{1} < x_{2}, then we also omit those a for which c^{2}x_{1} < a < c^{2}x_{2}.

Each point P found by this search is passed to the second procedure, which tests whether it has infinite order, discarding it if not. At the general stage we will have k independent points P_{i} for 1 ≤ i ≤ k (initially k = 0) which generate a subgroup A of rank k, and will have stored the k × k height pairing matrix M = (h(P_{i}, P_{j})) and its determinant R. Now we set P_{k+1} = P and compute ĥ(P_{i}, P_{k+1}) for i ≤ k + 1 to obtain a new height pairing matrix of order k + 1. If the determinant of this new matrix is non-zero, the new point is independent of the previous ones and we add it to the current list of generators, increment k, replace R by the new determinant, and go on with the point search. If the new determinant is zero, however, we use the values h(P_{i}, P) to express P_{k+1} as a linear combination of the P_{i} for i ≤ k, with approximate real coefficients: in fact we have

a_{1}P_{1} + a_{2}P_{2} + … + a_{k}P_{k} + a_{k+1}P_{k+1} = 0 (modulo torsion)

where for 1 ≤ i ≤ k + 1 the coefficient a_{i} is the (i, k + 1) cofactor of the enlarged matrix, which we will have stored during the computation of the new determinant. In particular, a_{k+1} is (up to sign) the previous value of R, and hence is non-zero. Next we find rational approximations to these floating-point coefficients a_{i} (using continued fractions, or MLLL if available), and clear denominators to obtain a new equation of the same form with coprime integer coefficients a_{i}, which we can check holds exactly. In this relation we still have a_{k+1} ≠ 0 (the first k points are independent). The simplest case now is when a_{k+1} = ±1, for then P_{k+1} is redundant and can be discarded. Similarly, if a_{i} = ±1 for some i ≤ k, then we may discard P_{i}, replacing it by P_{k+1}, and gaining index |a_{k+1}|. In general let a_{i} be the minimal non-zero coefficient (in absolute value); if |a_{i}| > 1, we find a coefficient a_{j} not divisible by a_{i} (which must exist since the coefficients are coprime) and write a_{j} = a_{i}q + b where 0 < b < |a_{i}|. Now since

a_{i}P_{i} + a_{j}P_{j} = a_{i}P_{i} + (a_{i}q + b)P_{j} = a_{i}(P_{i} + qP_{j}) + bP_{j}

we may replace the generator P_{i} by P_{i} + qP_{j}, replace the coefficient a_{j} by b (which is smaller than |a_{i}|), and replace i by j. After a finite number of steps we obtain a minimal coefficient a_{i} = 1 and can discard the current generator P_{i}, leaving a new set of k independent generators which generate a group larger than before by a finite index equal to the original value of |a_{k+1}|.

In this way, we will be able to find a ℤ-basis for the subgroup A of the Mordell–Weil group (modulo torsion) which is generated by the points of naive height less than the bound B. Often we know the rank r of our curve in advance, so that we can increase B until A has rank r. Then A has finite index in E(ℚ), and we must enlarge it to give the whole of E(ℚ). There are various methods one can use here, all of which rely on having explicit bounds for the difference between the naive and canonical heights on the curve E. The simplest general bound here is a result of Silverman (see [60]). One can certainly often obtain better bounds for individual curves, and there are also more complicated results which apply in general and which usually give much better bounds, such as the main result of [57].

For simplicity we will only give Silverman’s version of the bound. In the following proposition,^{3} the height of a rational number a/b with gcd(a, b) = 1 is h(a/b) = log max{|a|, |b|}, and log^{+}(x) = log max{1, |x|} for x ∈ ℝ.

Proposition 3.5.1. Let E be an elliptic curve defined by a standard Weierstrass equation over ℤ, with discriminant Δ and j-invariant j. Set 2^{∗} = 2 if b_{2} ≠ 0, or 2^{∗} = 1 if b_{2} = 0. Define

μ(E) = \frac{1}{6} (log |Δ| + log^{+}(j)) + log^{+}(b_{2}/12) + log(2^{∗}).

Then for all P ∈ E(ℚ),

−\frac{1}{12}h(j) − μ(E) − 1.922 ≤ ĥ(P) − h(P) ≤ μ(E) + 2.14.

This result is easiest to apply in the rank 1 case, as follows. Suppose we have a rational point P of infinite order on E, of height ĥ(P). If P is not a generator it is a multiple P = kQ (modulo torsion) of some generator Q, where k ≥ 2, so that ĥ(Q) ≤ \frac{1}{4} ĥ(P). By the preceding proposition we can bound the naive height of Q and adjust the bound B in our search accordingly. If a further search up to this bound finds no more points, then P was a generator after all; otherwise we are sure to find a generator.

Similar techniques are possible in higher rank situations, using estimates from the geometry of numbers. See the papers [70] and [57] for more details.

We may also remark that since P has finite order if and only if ĥ(P) = 0, the proposition implies that all torsion points have naive height h(P) ≤ \frac{1}{12}h(j) + μ(E) + 1.922, giving us another way of finding all the rational torsion points.

For the general case, the following simple result^{4} may be used.

Lemma 3.5.2. Let B > 0 be such that

S = {P ∈ E(ℚ) | ĥ(P) ≤ B}

contains a complete set of coset representatives for 2E(ℚ) in E(ℚ). Then S generates E(ℚ).

Proof. Let A be the subgroup of E(ℚ) modulo torsion generated by the points in S. Suppose that A is a proper subgroup; then we may choose Q ∈ E(ℚ) − A with ĥ(Q) minimal, since ĥ takes a discrete set of values. By hypothesis, there exist P ∈ A and R such that Q = P + 2R; certainly R ∉ A, so that ĥ(R) ≥ ĥ(Q) by minimality. Now using the fact that ĥ is quadratic and non-negative we obtain a contradiction:

ĥ(P) = \frac{1}{2}(ĥ(Q + P) + ĥ(Q − P)) − ĥ(Q)

≥ \frac{1}{2}ĥ(2R) − ĥ(Q)

= 2ĥ(R) − ĥ(Q) ≥ ĥ(Q) > B. □

We have two ways of using this in practice. First of all, it is possible to obtain from the two-descent procedure which we use to determine the rank (see the next section), a set of coset representatives for E(ℚ) modulo 2E(ℚ). Computing the heights of these points we can find a B for which the Lemma holds, to which we add the maximum difference between naive and canonical heights from the preceding proposition to get a bound on the naive heights of a set of generators.

3When referring to [60], recall that our ĥ is double Silverman’s; also, the constant 1.922 appearing here is a (normalized) correction, due to Bremner, of the constant in Silverman’s paper.

4Attributed in [60] to Zagier, it is also exercise 5 on page 84 of Cassels’ book [8].

Alternatively, assuming that we know the rank r, we first run our search until we find r independent points P_{i}. Now it is easy to check whether a point P is twice another: if any subset of the P_{i} sums to 2Q for some Q we replace one of the P_{i} in the sum by Q and gain index 2. After a finite number of steps (since we are in a finitely-generated group) we obtain independent points which are independent modulo 2, and proceed as before.

Again, we have only presented here the most straightforward strategies for enlarging a set of r independent points in E(ℚ) to a full ℤ-basis; this is a topic of active research, with new ideas being developed rapidly: see the paper [57] for some recent advances.

Putting the pieces together, we can determine a set of generators for E(ℚ) modulo torsion, and then compute the regulator, provided that we know its rank. If we do not know the rank, we at least can obtain lower bounds for the rank. Together with the torsion points found in section 3.3, we will have determined the Mordell–Weil group E(ℚ) explicitly. Computing the rank is the subject of the next section.

3.6 The Mordell–Weil group III: the rank

For an elliptic curve E defined over the rationals, the rank of the Mordell–Weil group E(ℚ) is by far the hardest of the elementary quantities associated with E to compute, both theoretically and in terms of implementation. Strictly speaking, the two-descent algorithms we will describe are not algorithms at all, as they are not guaranteed to terminate in all cases. One part of the procedure involves establishing whether or not certain curves of genus one have rational points, when they are known to have points everywhere locally (that is, over ℝ and over the p-adic field ℚ_{p} for all primes p): there is no known algorithm to decide this in general. Moreover, even without this difficulty, for curves with large coefficients and no rational points of order two, the general two-descent algorithm takes too long to run in practice. For simplicity, we will refer to the procedures as rank algorithms, although their output in certain cases will be bounds on the rank rather than its actual value.

We originally decided to implement a general two-descent procedure in order to check that the modular curves we had computed did have their rank equal to the analytic rank, which we knew, as described in the previous chapter. This was a somewhat thankless task, as it involved a large programming effort, and a large amount of computer time to run the resulting program, in order to verify that approximately 2500 numbers did in fact have the values 0, 1 or 2 which we were already sure were correct. Since the project started, the major theoretical advances by Kolyvagin, Rubin and others meant that all the cases of rank 0 or 1 were known anyway, which left just 18 cases of conjectured rank 2 to verify. In the end we were able to verify these cases, and to check all but a few dozen of the rank 0 or 1 curves; we also obtained extra information by the two-descent procedure, such as the 2-rank of the Tate–Shafarevich group � , and a set of coset representatives for E(ℚ)/2E(ℚ).

Since the original implementation, the algorithm has been much improved in many ways (notably the syzygy sieve in the search for quartics, the systematic use of group structure in the 2-isogeny case, and the use of quadratic sieving in searching for rational points on homogeneous spaces: see below for details). Our program mwrank,^{5} based on the algorithm, now works well on a much larger set of curves, including some of fairly high rank such as a curve of Fermigier [23] with rank 13 and 2-torsion (see the example below), and several curves with no 2-torsion and ranks 6, 7 and 8. However, curves with extremely large coefficients, such as Nagao’s curve of rank (at least) 21 (see [45]), are beyond the reach of this algorithm owing to the enormous search regions required. One can also use the program mwrank to find points on curves which are too large to find by the search methods of the previous section.

5Available from the author’s ftp site: see the Introduction for details.

We will not describe here the theory of two-descent, which is the basis of the algorithm, in great detail. Roughly speaking, one has an injective homomorphism from E(ℚ)/2E(ℚ) into a finite elementary abelian 2-group, the 2-Selmer group, and attempts to determine the image; if this has order 2^{t} then the rank of E(ℚ) is t, t − 1 or t − 2 according to whether the number of points of order 2 in E(ℚ) is 0, 1 or 3 (respectively). This procedure applies to arbitrary curves, and is called general two-descent. When E has a rational point P of order two, there is a rational 2-isogeny ϕ : E → E′ = E/ ⟨P ⟩ and a dual isogeny ϕ′ : E′ → E. We may then proceed differently, using a procedure we call two-descent via 2-isogeny: we embed each of E′/ϕ(E) and E/ϕ′(E′) into finite subgroups of ℚ^{∗}/(ℚ^{∗})^{2}, which are easy to write down. This is in contrast to the general two-descent, where one has to work hard to find the Selmer group itself. A full description of two-descent can be found in the standard references such as the books by Silverman [58], Husemöller [27], Knapp [28], or Cassels [8], but the descriptions given there are only easy to apply when E has all its 2-torsion rational. For the general case where there are no rational points of order 2, the main reference is one of the original papers [3] by Birch and Swinnerton–Dyer on their Conjecture, and we followed that paper closely in writing the first version of our program. More detail on the invariant theory, which has resulted in substantial improvements to the general two-descent algorithm, can be found in the paper [20]; a very full description of the algorithm, together with its extension to real quadratic number fields (see also [19]), can be found in Serf’s thesis [52].

Both algorithms involve the classification of certain curves, associated with the given elliptic curve E, called principal homogeneous spaces. These are twists of E: curves of genus 1 isomorphic to E over an extension field, but not (necessarily) over ℚ itself; they need not have rational points, so need not themselves be elliptic curves. When they do have rational points, these map to rational points on E; the maps H → E are called 2-coverings and have degree 4 (in the general two-descent) or 2 (in the 2-isogeny descent). The homogeneous spaces which arise in both algorithms have equations of the form

(3.6.1)  H :  y^{2} = g(x) = ax^{4} + bx^{3} + cx^{2} + dx + e

where g(x) is a quartic polynomial with rational coefficients. For brevity we will usually refer to these principal homogeneous spaces simply as quartics. The invariants I and J of g(x) (see below for their definition) are related to the invariants c_{4} and c_{6} of either E or the 2-isogenous curve E′. In the case of descent via 2-isogeny, g(x) will in fact be a quadratic in x^{2}. We will be interested in whether the quartic H has points over ℚ or one of its completions, the p-adics ℚ_{p} or the reals ℝ. Such a point will either be an affine point (x, y) satisfying the equation (3.6.1), or one of the two points at infinity on the projective completion of H, which are rational if and only if a is a square.

In all cases, a quartic with a (global) rational point (x, y) will lead to a rational point on the original curve E, and the set of all the rational points thus obtained will cover the cosets of 2E(ℚ) in E(ℚ); thus we will be able to determine the rank of E(ℚ), and at the same time obtain a set of points which generates a subgroup of the Mordell-Weil group E(ℚ) of odd, finite index. Quartics with no global rational point which are everywhere locally soluble arise from non-trivial elements in the Tate–Shafarevich group of E (or of E′); if these exist, we will only obtain upper and lower bounds for the rank. This is because we currently have no general procedure for proving that a quartic with no rational points does have none. In practice, moreover, it is often impossible to distinguish between such a quartic and one with rational points which are all very large, and hence outside the search region; this happens when a curve has some very large generators, and in such cases also we may only be able to give bounds for the rank. Further work on these questions is clearly needed, and is currently the focus of much active research.

Since the covering maps H → E have degree 2 or 4, the rational points on H tend to be smaller (in the sense of naive height) than the rational points they map to on E; this makes them easier to find by search. Here is an example of this: the curve y^{2} = x^{3} − 673 has rank 2, with generators P_{1} = (29, 154) and P_{2} = (33989323537/61761^{2}, −1384230292401340/61761^{3}). The second generator, which would take a very long time to find by searching on the curve itself, is obtained from the rational point (x, y) = (191/97, 123522/97^{2}) on the quartic with coefficients (a, b, c, d, e) = (−2, 4, −24, 164, −58). This is much easier to find: our program takes less than a second to find the rank and both generators of this curve (but in this time it does not prove that they are generators, only that they generate a subgroup of finite odd index in the Mordell-Weil group).

Before we describe the two main two-descent algorithms, we will present algorithms for determining local solubility and for attempting to determine global solubility of a quartic equation such as (3.6.1), as these are used in both the algorithms.

Checking local solubility.

Here we present an algorithm for determining the local solubility of a curve of the form (3.6.1), where g(x) is a square-free quartic polynomial with integer coefficients. It is easy to generalize this algorithm in two ways: firstly, one might be interested in polynomials of higher degree (when studying curves of higher genus, for example); secondly, working over a general number field K, one would replace the p-adic field ℚ_{p} here with the appropriate completion of K. These extensions are quite straightforward.

Solubility at the infinite prime (that is, over the reals) is easily determined. If g(x) has a real root then it certainly takes positive values, so that H has real points; if g(x) has no real roots, then the values of g(x) have constant sign, and we merely have to check that a > 0.

Regarding the finite primes, we first observe that there are only a finite number which need checking in each case, for if p is an odd prime not dividing the discriminant of g, then H certainly has points modulo p which are nonsingular and hence lift to p-adic points. For the other primes, we present an algorithm first given in [3].

It suffices to determine solubility in ℤ_{p} for either g(x) or g^{∗}(x) = ex^{4} + dx^{3} + cx^{2} + bx + a, and in the latter case we may assume x ∈ pℤ_{p}. Given x_{k} modulo p^{k}, one tries to lift to a p-adic point (x, y) with x ≡ x_{k} (mod p^{k}). In [3], conditions are given for this to be possible; more precisely, one of three possibilities may occur (given k and x_{k}): either a lifting is definitely possible, and we may terminate the algorithm with a positive result; or it is definitely not possible, and we reject this value of x_{k}; or it is impossible to decide without considering x_{k} modulo a higher power of p. The test for this lifting is given below in the two subroutines called lemma6 and lemma7, named after the corresponding results in [3]. This leads to a recursive algorithm which is guaranteed to terminate since in any given case there is an exponent k such that it is possible^{6} to determine p-adic solubility by considering solubility modulo p^{k}. All this is an exercise in Hensel’s Lemma; the prime p = 2 needs to be considered separately. For the details, we refer to the pseudocode below, or to [3]. Further information on local solubility may be found in [56] and [57].

Here is the pseudocode for these algorithms. Note that for any given elliptic curve, all the homogeneous spaces considered will have the same discriminant as the curve (up to a power of 2), so that in practice we would not need to factorize the discriminant of each quartic.

6In fact, if k > ord_{p}(disc(g)), and also k ≥ 2 when p = 2, then the third possibility cannot occur in algorithms lemma6 and lemma7.

Algorithm for determining local solubility of a quartic

INPUT:        a, b, c, d, e (integer coefficients of a quartic g(x))
OUTPUT:       TRUE/FALSE    (solubility of y2=g(x) in ℝ and in ℚ_{p} for all p)
1.  BEGIN
2.  IF NOT R soluble(a,b,c,d,e) THEN RETURN FALSE FI;
3.  IF NOT Qp soluble(a,b,c,d,e,2) THEN RETURN FALSE FI;
4. Δ = discriminant(a,b,c,d,e);
5.  p list = odd prime factors(Δ);
6.  FOR p IN p list DO
7.  BEGIN
8.       IF NOT Qp soluble(a,b,c,d,e,p) THEN RETURN FALSE FI
9.  END;
10.  RETURN TRUE
11.  END

(Subroutine for determining real solubility)

SUBROUTINE R soluble(a,b,c,d,e)

INPUT:        a, b, c, d, e (integer coefficients of a quartic g(x))
OUTPUT:       TRUE/FALSE    (solubility of y2=g(x) in ℝ)
1.  BEGIN
2.  IF a>0 THEN RETURN TRUE FI;
3.  x list = real roots(a*x4+b*x3+c*x2+d*x+e=0);
4.  IF length(x list)>0 THEN RETURN TRUE FI;
5.  RETURN FALSE
6.  END

(Subroutine for determining p-adic solubility)

SUBROUTINE Qp soluble(a,b,c,d,e,p)

INPUT:        a, b, c, d, e (integer coefficients of a quartic g(x))
             p             (a prime)
OUTPUT:       TRUE/FALSE    (solubility of y2=g(x) in ℚ_{p})
1.  BEGIN
2.  IF Zp soluble(a,b,c,d,e,0,p,0) THEN RETURN TRUE FI;
3.  IF Zp soluble(e,d,c,b,a,0,p,1) THEN RETURN TRUE FI;
4.  RETURN FALSE
5.  END

(Recursive ℤ_{p}-solubility subroutine)

SUBROUTINE Zp soluble(a,b,c,d,e,x k,p,k)

INPUT:        a, b, c, d, e (integer coefficients of a quartic g(x))
             p             (a prime)
             x k           (an integer)
             k             (a non-negative integer)
OUTPUT:       TRUE/FALSE    (solubility of y2=g(x) in ℤ_{p}, with x≡x k (mod p^{k}))
1.  BEGIN
2.  IF p=2
3.  THEN code = lemma7(a,b,c,d,e,x k,k)
4.  ELSE code = lemma6(a,b,c,d,e,x k,p,k)

5.  FI;
6.  IF code=+1 THEN RETURN TRUE FI;
7.  IF code=-1 THEN RETURN FALSE FI;
8.  FOR t = 0 TO p-1 DO
9.  BEGIN
10.      IF Zp soluble(a,b,c,d,e,x k+t*p^{k},p,k+1) THEN RETURN TRUE FI
11.  END;
12.  RETURN FALSE
13.  END

(ℤ_{p} lifting subroutine: odd p)

SUBROUTINE lemma6(a,b,c,d,e,x,p,n)

1.  BEGIN
2.  gx = a*x4+b*x3+c*x2+d*x+e;
3.  IF p adic square(gx,p) THEN RETURN +1 FI;
4.  gdx = 4*a*x3+3*b*x2+2*c*x+d;
5.  l = ord(p,gx); m = ord(p,gdx);
6.  IF (l≥m+n) AND (n>m) THEN RETURN +1 FI;
7.  IF (l≥2*n) AND (m≥n) THEN RETURN 0 FI;
8.  RETURN -1
9.  END

(ℤ_{2} lifting subroutine)

SUBROUTINE lemma7(a,b,c,d,e,x,n)

1.  BEGIN
2.  gx = a*x4+b*x3+c*x2+d*x+e;
3.  IF p adic square(gx,2) THEN RETURN +1 FI;
4.  gdx = 4*a*x3+3*b*x2+2*c*x+d;
5.  l = ord(p,gx); m = ord(p,gdx);
6.  gxodd = gx; WHILE even(gxodd) DO gxodd = gxodd/2;
7.  gxodd = gxodd (mod 4);
8.  IF (l≥m+n) AND (n>m) THEN RETURN +1 FI;
9.  IF (n>m) AND (l=m+n-1) AND even(l) THEN RETURN +1 FI;
10.  IF (n>m) AND (l=m+n-2) AND (gxodd=1) AND even(l) THEN RETURN +1 FI;
11.  IF (m≥n) AND (l≥2*n) THEN RETURN 0 FI;
12.  IF (m≥n) AND (l=2*n-2) AND (gxodd=1) THEN RETURN 0 FI;
13.  RETURN -1
14.  END

A few further remarks on these algorithms: firstly, only trivial changes need to be made to the algorithms Qp soluble and Zp soluble to make them apply to more general equations of the form y^{2} = g(x) where g(x) is a non-constant squarefree integer polynomial. This is relevant for work on curves of higher genus, and was observed by S. Siksek. Secondly, extensions to more general p-adic fields are also useful in studying curves over number fields, and again the extensions of Lemma 6 and Lemma 7 in [3] are not difficult. See the theses [56] and [52] for details of such extensions.

Lastly, D. Simon observed that in our application of the algorithms lemma6 and lemma7, we only care whether there is a solution, not necessarily that there is a solution congruent to the given x (mod p^{k}); hence line 6 of subroutine lemma6 and line 8 of subroutine lemma7 can both be replaced by:

    IF l>2*m THEN RETURN +1 FI.

Checking global solubility.

To determine whether an equation (3.6.1) has a rational point is much harder than the corresponding local question. All we can do at present is search (efficiently) for a point up to a certain height, after checking that there is no local obstruction. The only satisfactory way known at present to decide on the existence of rational points on these homogeneous spaces is to carry out so-called higher descents; as mentioned above, this is the subject of current work (see [63], for example), and we will not consider it further here.

Our strategy is to look first for a small rational point, using a very simple procedure with low overheads; if this fails, we check for local solubility; if this passes, we start a much more thorough search for a global point, using a quadratic sieving procedure rather similar to the one described in the previous section for finding points on the elliptic curve itself. (In fact, such a sieve-assisted search may be used to find rational points on any curve given by an equation of the form y^{2} = g(x) where g(x) is a polynomial in x.) The philosophy here is that there is no point in looking hard for rational points unless one is sure of local solubility, but also that there is no point in checking local solubility when there is an obvious global point.

To carry out the sieve-assisted search, for each possible denominator of x one precomputes, for each of several sieving moduli m, the residues to which the numerator of x must belong if the right-hand side of the equation is to be a square modulo m. In addition, it is easy to see that for every odd prime p dividing the denominator of the x-coordinate of a rational point, we must have (\frac{a}{p}) = +1; so provided that the leading coefficient a is not a square (in which case the points at infinity are rational anyway), we precompute a list of primes p for which (\frac{a}{p}) = −1, and discard possible denominators divisible by any of these primes. For p = 2 a similar condition holds.^{7} One also obviously restricts the search to ranges of x for which g(x) is positive; depending on the number of real roots of g and the sign of a, this splits the search into up to three intervals. Finally, in the case of two-descent via 2-isogeny, where the quartics are polynomials in x^{2} and thus even, we may restrict to positive x.

For reasons of space, we will only give here the code for a simple point search with no sieving.

Algorithm for searching for a rational point on a quartic: simple version

INPUT:        a, b, c, d, e (integer coefficients of a quartic g(x))
              k1, k2        (lower and upper bounds)
OUTPUT:       TRUE/FALSE    (solubility of y2=g(x) in ℚ with x=u/w
                            and k1 ≤ |u|+w ≤ k2)
 1.  BEGIN
 2.  FOR n = k1 TO k2 DO
 3.  BEGIN
 4.       IF n=1 THEN
 5.           IF square(a) RETURN TRUE FI;
 6.           IF square(e) RETURN TRUE FI
 7.       ELSE
 8.           FOR u = 1 TO n-1 DO
 9.           BEGIN
10.               IF gcd(u,n)=1
11.               THEN
12.                    w = n-u;
13.                    IF square(a*u4+b*u3w+c*u2w2+d*uw3+e*w4) RETURN TRUE FI;

7I am grateful to J. Gebel for this idea, which saves considerable time in practice.

14.                    IF square(a*u4-b*u3w+c*u2w2-d*uw3+e*w4) RETURN TRUE FI
15.               FI
16.           END
17.      FI
18.  END;
19.  RETURN FALSE
20.  END

We will now describe the two main two-descent algorithms: two-descent via 2-isogeny for use when E has a rational point of order 2, and general two-descent in the general case. We only use general two-descent when there is no point of order 2, so that the first method does not apply. The situation is not appreciably simpler when E has all three of its points of order two rational than when there is just one rational point of order two, and so we will not bother to consider this case separately.

Method 1: descent using 2-isogeny.

Suppose that E has a rational point P of order 2. By a change of coordinates we may assume that E has equation

E : y^{2} = x(x^{2} + cx + d)

where P = (0, 0), and c, d ∈ ℤ. Explicitly, in terms of a Weierstrass equation, let x_{0} be a root of the cubic x^{3} + b_{2}x^{2} + 8b_{4}x + 16b_{6}, and set c = 3x_{0} + b_{2}, d = (c + b_{2})x_{0} + 8b_{4}. If a_{1} = a_{3} = 0, then we can avoid a scaling factor of 2 by letting x_{0} be a root of x^{3} + a_{2}x^{2} + a_{4}x + a_{6}, and setting c = 3x_{0} + a_{2}, d = (c + a_{2})x_{0} + a_{4}. The 2-isogenous curve E′ = E/ ⟨P ⟩ has equation

E′ : y^{2} = x(x^{2} + c′x + d′)

where

c′ = −2c  and  d′ = c^{2} − 4d.

The nonsingularity condition on E is equivalent to dd′ ≠ 0. The 2-isogeny ϕ: E → E′ has kernel \{0, P \} and in general maps (x, y) to (\frac{y^{2}}{x^{2}} , \frac{y(x^{2}−d)}{x^{2}}). The dual isogeny ϕ′: E′ → E maps

(x, y) to (\frac{y^{2}}{4x^{2}} , \frac{y(x^{2}−d′)}{8x^{2}}).

For each factorization d = d_{1}d_{2}, with d_{1} square-free, we consider the homogeneous space

H(d_{1}, c, d_{2}) : v^{2} = d_{1}u^{4} + cu^{2} + d_{2}.

Let n_{1} = n_{1}(c, d) be the number of factorizations of d for which the quartic H(d_{1}, c, d_{2}) has a rational point, and n_{2} = n_{2}(c, d) the number for which the quartic has a point everywhere locally. Define n′_{1} = n_{1}(c′, d′) and n′_{2} = n_{2}(c′, d′) similarly. Then it is not hard to show by rather explicit calculation (see below and the references given) that E(ℚ)/ϕ′(E′(ℚ)) is isomorphic to the subgroup of ℚ^{∗}/(ℚ^{∗})^{2} generated by the factors d_{1} for which H(d_{1}, c, d_{2}) has a rational point. Thus

|E(ℚ)/ϕ′(E′(ℚ))| = n_{1},

which must therefore be a power of 2, say n_{1} = 2^{e_{1}}; similarly,

|E′(ℚ)/ϕ(E(ℚ))| = n′_{1} = 2^{e′_{1}}.

It then follows (see below) that

(3.6.2)  rank(E(ℚ)) = rank(E′(ℚ)) = e_{1} + e′_{1} − 2.

With luck one will find rational points on all the quartics which have them everywhere locally; then n_{1} = n_{2}, and there is no ambiguity in the result. However there will be cases in which the number ñ_{1} of quartics on which we can find a rational point is strictly less than n_{2}. In such cases, we will only have upper and lower bounds for n_{1}, and similarly for n′_{1}, leading to upper and lower bounds for the rank. This can happen for two reasons: either there is a rational point on some quartic, but our search bound was too small to find it; or the quartic has points everywhere locally but no global rational point.

The quartics H which have points everywhere locally but not globally come from elements of order 2 in the Tate–Shafarevich groups � (E/ℚ) and � (E′/ℚ). There is an exact sequence

0 → E(ℚ)/ϕ′(E′(ℚ)) → S^{(ϕ′)}(E′/ℚ) → � (E′/ℚ)[ϕ′] → 0

coming from Galois cohomology; here S^{(ϕ′)}(E′/ℚ) is the Selmer group of order n_{2} whose elements are represented by the homogeneous spaces H(d_{1}, c, d_{2}) which are everywhere locally soluble, and � (E′/ℚ) is the Tate–Shafarevich group of E′. The injective map E(ℚ)/ϕ′(E′(ℚ)) → S^{(ϕ′)}(E′/ℚ) is induced by taking a point (x, y) in E(ℚ) with x ≠ 0 to the space H(d_{1}, c, d_{2}) where d_{1} = x modulo squares: if x = d_{1}u^{2} and v = uy/x then (u, v) is a rational point on H(d_{1}, c, d_{2}). The point P = (0, 0) maps to d modulo squares. Conversely, if (u, v) is a rational point on H(d_{1}, c, d_{2}) then (x, y) = (d_{1}u^{2}, d_{1}uv) is a rational point on E. (In proving these statements, one has to check that two rational points on E have the same x-coordinate modulo squares if and only if their difference is in ϕ′(E′(ℚ)); for example, the image of P is d, which is a square if and only if P ∈ ϕ′(E′(ℚ)).) It follows that n_{1} is the order of E(ℚ)/ϕ′(E′(ℚ)), as stated above, and hence that

|� (E′/ℚ)[ϕ′]| = n_{2}/n_{1}.

Similarly, from the exact sequence

0 → E′(ℚ)/ϕ(E(ℚ)) → S^{(ϕ)}(E/ℚ) → � (E/ℚ)[ϕ] → 0

with similarly defined maps, we obtain

|� (E/ℚ)[ϕ]| = n′_{2}/n′_{1}.

Thus the result is only genuinely ambiguous when either � (E/ℚ)[ϕ] or � (E′/ℚ)[ϕ′] is non-trivial, so that not all elements of the Selmer groups are obtained from rational points on the elliptic curves. This is rare for the curves in the tables, but obviously must be taken into account in general. A typical situation is to have n_{2}n′_{2} = 16 and n_{1}n′_{1} ≥ 4, when one suspects that r = 0 with |� (E/ℚ)[2]| = 4 or |� (E′/ℚ)[2]| = 4, but where it is possible instead that r = 2 and |� (E/ℚ)[2]| = |� (E′/ℚ)[2]| = 1. Curve 960D1 in the tables is an example of this, although in this case since the curve is modular and we know that L(E, 1) ≠ 0, it must have rank 0 by the result of Kolyvagin mentioned earlier. We can also deduce this by working with the 2-isogenous curves 960D3 and 960D2, where there is no ambiguity: here n_{1} = n_{2} = n′_{1} = n′_{2} = 2, showing that the rank is certainly 0. (Note that isogenous curves have the same rank, but not necessarily the same order of � , which can work to our advantage in cases like this.) Returning to the pair 960D1-960D2 where we compute n_{1} = n_{2} = 1, n′_{2} = 16 and n′_{1} ≥ 4, now we know that the rank is in fact zero we can conclude that n′_{1} = 4, and that |� (E/ℚ)[ϕ]| = 4. The nontriviality of � (E/ℚ) in this case is confirmed by the Birch–Swinnerton-Dyer conjecture, which for this curve predicts that � has order 4 (see Table 4).

Local solubility of H(d_{1}, c, d_{2}) is automatic for all primes p which do not divide 2dd′; for those p which do divide 2dd′ we may apply the general criteria of Birch and Swinnerton-Dyer.

Local solubility in ℝ is easy to determine here: if d′ < 0 then we require d_{1} > 0, while if d′ > 0 then either d_{1} > 0 or c + \sqrt{d′} > 0 is necessary. Thus if either d′ < 0, or d′ > 0 and c + \sqrt{d′} < 0, then we only consider positive divisors d_{1} of d, and need not apply the general test for solubility in ℝ.

Each rational point (u, v) on H(d_{1}, c, d_{2}) maps, as observed above, to the point (d_{1}u^{2}, d_{1}uv) on E; modulo ϕ′(E′(ℚ)), this is independent of the rational point (u, v), and only depends on d_{1} modulo squares. Similarly, a rational point (u, v) on H(d′_{1}, c′, d′_{2}) maps to a point on E′, and hence via the dual isogeny ϕ′ to the point

(\frac{v^{2}}{4u^{2}} , \frac{v(d′_{1}u^{4} − d′_{2})}{8u^{3}})

in E(ℚ). The set of n_{1}n′_{1} points in E(ℚ) thus determined (by adding the points constructed in this way) cover the cosets of E(ℚ)/2E(ℚ), either once each, when |E(ℚ)[2]| = 4, which is when d′ is a square, or twice, when d′ is not a square and |E(ℚ)[2]| = 2. Thus, when |E(ℚ)[2]| = 2 we have

\frac{n_{1}n′_{1}}{2} = |E(ℚ)/2E(ℚ)| = 2^{r+1},

while if |E(ℚ)[2]| = 4 we have

n_{1}n′_{1} = |E(ℚ)/2E(ℚ)| = 2^{r+2};

hence 2^{r} = n_{1}n′_{1}/4 in both cases, proving (3.6.2).

When counting n_{1} and n_{2} (and similarly, n′_{1} and n′_{2}), it is very useful to use the fact that each is a power of 2, being the order of an elementary abelian 2-group. This is particularly important when d (or d′) has many distinct prime factors. Let A_{0} be the group of all divisors of d modulo squares, of order n_{0} (say). Then A_{0} is generated by −1 and the primes dividing d, so that n_{0} = 2^{e_{0}} where e_{0} is the number of distinct prime factors of d, plus 1. Within A_{0} we must determine the subgroups A_{1} and A_{2} of orders n_{1} and n_{2}, consisting of those divisors d_{1} of d for which the corresponding homogeneous space is everywhere locally or globally soluble, respectively.

We can effectively reduce the size of the set A_{0} of divisors to be searched by a factor up to 8 as follows: as observed above, if either d′ < 0, or d′ > 0 and c + \sqrt{d′} < 0, then we need only consider positive divisors d_{1} of d, cutting in two the number of elements of A_{0} which may lie in A_{1}. Secondly, we may take advantage of the fact that we know the rational point (0, 0) on E; thus we know that d is in A_{2} (though possibly just the identity if d is a square); similarly, if d′ is a square then x^{2} + cx + d factorizes, say as (x − x_{2})(x − x_{3}), and we know that x_{2} and x_{3} also lie in A_{2}.

More generally, whenever we find in the course of our systematic search through the elements of A_{0} that the element d_{1} lies in A_{2}, we can effectively factor out d_{1} and reduce the number of remaining values to check by a factor of 2. Of course, this requires careful book-keeping in the implementation; for simplicity, we omit these refinements from the pseudocode below, where we simply loop over all square-free divisors of d and d′.

As an example of the saving that can be made, consider the curve of rank 13 constructed by Fermigier in [23]; this is of the form y^{2} = x(x^{2} + cx + d) with

c = 36861504658225  and
d = 1807580157674409809510400 = 2^{15} · 3^{4} · 5^{2} · 7^{2} · 17 · 23 · 29 · 41 · 103 · 113 · 127 · 809,

so that d has 12 distinct prime factors and 2^{13} = 8192 square-free divisors. Since d is non-square we can cut the set in half, say by excluding all d_{1} divisible by the largest prime factor 809, leaving 4096 values to test. In our implementation, the results of the test are as follows:

• 7 non-trivial values of d_{1} give rational points after searching, as well as d_{1} = 1 which gives the trivial point;

• 120 further values are in the subgroup A_{2} generated by these 7 values and need not be tested;

• 122 further values were tested and found to be not everywhere locally soluble, hence not in A_{1};

• 3846 further values were discarded as being a product of an element of A_{2} and an element not in A_{1}, and hence not in A_{1}.

Thus in this case we find that n_{1} = n_{2} = 256, after only having to search for points on seven homogeneous spaces. Working with the isogenous curve, we obtain n′_{1} = n′_{2} = 128 after only searching six homogeneous spaces for points. Thus e_{1} = 8, e′_{1} = 7 and the rank is 13. Note that in the course of computing this value, we have searched precisely 13 homogeneous spaces, and the points we thereby construct give 13 generators of E(ℚ)/2E(ℚ) modulo torsion. Adding P = (0, 0) to this list gives 14 points which generate E(ℚ)/2E(ℚ) (which has order 2^{14}), and which therefore generate a subgroup of finite odd index in the full Mordell-Weil group E(ℚ).

The situation is not always this simple, however, even for curves where � [2] is trivial, since there may be homogeneous spaces with rational points which are hard to find. For example, consider Fermigier’s curve of rank 14 from [23], with c = 2429469980725060 and d = 275130703388172136833647756388 (which has 14 prime factors). When we run our program using a (logarithmic) bound of 10 in the search for rational points on the quartics, we find n_{1} ≥ 64, n′_{1} ≥ 128, while n_{2} = n′_{2} = 256. Here the correct values are n_{1} = n′_{1} = 256, giving r = 14, but we only find 11 ≤ r ≤ 14; and in the process, we have had to search many more homogeneous spaces for rational points.

Here is the pseudo-code which implements the algorithm just described. The main routine aborts if either the input curve is singular (this is useful if one wants to apply the algorithm systematically to a range of inputs) or if there is no point of order two. The latter is detected in lines 6–7, where an integer root to a monic cubic with integer coefficients is found (if it exists). Most of the work is done in the subroutine count(c,d,p list) which determines n_{2}(c, d) and, as far as possible, n_{1}(c, d). Here p list is the set of ‘bad’ primes dividing 2dd′ where local solubility needs to be checked, which we only compute once. There are two calls to the subroutine rational point(a,b,c,d,e,k1,k2), which seeks a rational u/w with k_{1} ≤ |u| + w ≤ k_{2} such that g(u/w) is a rational square, where g(x) = ax^{4} + bx^{3} + cx^{2} + dx + e. (Here w > 0 and gcd(u, w) = 1.) In the first call we carry out a quick check for ‘small’ points; then we look further, having first checked for everywhere local solubility. The particular parameters lim1, lim2 for the search will probably be decided at run time. The subroutines Qp soluble and rational point are implementations of the algorithms given earlier (though in practice we would use a more efficient algorithm for the second call to rational point, as explained above).

Algorithm for computing rank: rational 2-torsion case

INPUT:        a1, a2, a3, a4, a6(coefficients of E)
OUTPUT:       r min, r max       (bounds for rank of E)
              S,S'              (upper bounds for #� (E)[ϕ] and #� (E′)[ϕ′])
 1.  BEGIN
 2.  IF a1=a3=0
 3.  THEN s2 = a2; s4 = a4; s6 = a6
 4.  ELSE s2 = a1*a1+4*a2; s4 = 8*(a1*a3+2*a4); s6 = 16*(a3*a3+4*a6)
 5.  FI;
 6.  x list = integer roots(x3+s2*x2+s4*x+s6=0);

 7.  IF length(x list)=0 THEN abort ELSE x0 = x list[1] FI;
 8.  c = 3*x0+s2; d = (c+s2)*x0 + s4;
 9.  c' = -2*c; d' = c2-4*d;
10.  IF d*d'=0 THEN abort FI;
11.  p list = prime divisors(2*d*d');
12.  (n1,n2) = count(c,d,p list);
13.  (n1',n2') = count(c',d',p list);
14.  e1 = log 2(n1); e2 = log 2(n2);
15.  e1' = log 2(n1'); e2' = log 2(n2');
16.  r min = e1+e1'-2; r max = e2+e2'-2;
17.  S = n2'/n1'; S' = n2/n1;
18.  RETURN r min, r max, S, S'
19.  END

(Main counting subroutine)

SUBROUTINE count(c,d,p list)

 1.  BEGIN
 2.  n1 = n2 = 1; d' = c2-4*d;
 3.  d1 list = squarefree divisors(d);
 4.  FOR d1 IN d1 list DO
 5.  BEGIN
 6.       IF rational point(d1,0,c,0,d/d1,1,lim1)
 7.       THEN n1 = n1+1; n2 = n2+1
 8.       ELSE
 9.           IF everywhere locally soluble(c,d,d',d1,p list)
10.           THEN
11.               n2 = n2+1;
12.               IF rational point(d1,0,c,0,d/d1,lim1+1,lim2)
13.               THEN n1 = n1+1
14.               FI
15.           FI
16.      FI
17.  END;
18.  RETURN (n1, n2)
19.  END

(Subroutine to check for everywhere local solubility)

 1.  SUBROUTINE everywhere locally soluble(c,d,d',d1,p list)
 2.  BEGIN
 3.  IF d'<0 AND d1<0 THEN RETURN FALSE FI;
 4.  IF d'>0 AND d1<0 AND (c+sqrt(d'))<0 THEN RETURN FALSE FI;
 5.  FOR p IN p list DO
 6.  BEGIN
 7.       IF NOT Qp soluble(d1,0,c,0,d/d1,p) THEN RETURN FALSE FI
 8.  END;
 9.  RETURN TRUE
10.  END

Method 2: general two-descent.

We now turn to the general two-descent, which applies whether or not E has a rational point of order 2. Again, the basic idea is to associate to E a collection of 2-covering quartic curves (or homogeneous spaces) H. These have equations of the form

(3.6.1)  H :  y^{2} = g(x) = ax^{4} + bx^{3} + cx^{2} + dx + e

with a, b, c, d, e ∈ ℚ, such that the invariants

I = 12ae − 3bd + c^{2}  and  J = 72ace + 9bcd − 27ad^{2} − 27eb^{2} − 2c^{3}

are related to the c_{4} and c_{6} invariants of E via

I = λ^{4}c_{4}  and  J = 2λ^{6}c_{6}

for some λ ∈ ℚ^{∗}. Two such quartics g_{1}(x), g_{2}(x) are equivalent if

g_{2}(x) = μ^{2}(γx + δ)^{4}g_{1}(\frac{αx + β}{γx + δ} )

for some α, β, γ, δ and μ ∈ ℚ, with μ and αδ − βγ non-zero. The invariants of g_{1}(x) and g_{2}(x) are then related by the scaling factor λ = μ(αδ − βγ):

I(g_{2}) = μ^{4}(αδ − βγ)^{4}I(g_{1}),
J(g_{2}) = μ^{6}(αδ − βγ)^{6}J(g_{1}).

We set Δ = 4I^{3} − J^{2} = 27disc(g), and call Δ the discriminant.

In particular, by scaling up the coefficients, we may assume that the invariants I and J are integral. The number of equivalence classes of quartics with given invariants (up to a scaling factor λ) which are everywhere locally soluble is finite. One of our tasks will be to determine, for a given integral quartic, an equivalent integral one with minimal invariants. This process is closely analogous to the one considered earlier in this chapter, using Kraus’s conditions or Tate’s algorithm to determine minimal models for elliptic curves. Indeed, we will see below that if c_{4} and c_{6} are invariants of a minimal model for the elliptic curve E, then I = c_{4} and J = 2c_{6} are also minimal, except possibly at the prime 2. (We may lose minimality at 2 because the equations (3.6.1) we use for homogeneous spaces are not completely general, not having terms in y, xy or x^{2}y; to remove these by completing the square involves a scaling by a factor of 2.)

We now explain the relationship between equivalence classes of soluble quartics with invariants I and J and rational points on the elliptic curve. More details of this relationship, including proofs, may be found in [20]. For convenience, we again start by making a coordinate transformation: if c_{4} and c_{6} are the integral invariants of our curve E, we set I = c_{4} and J = 2c_{6}, and replace E by the isomorphic curve

(3.6.3)  E_{I,J} : Y^{2} = F (X) = X^{3} − 27IX − 27J.

This is the model on which the rational points we construct will naturally lie; it is then a simple matter to transfer them back to the original model for E. For simplicity, we will still continue to refer to the curve simply as E when this will not cause confusion.

Associated to each quartic g there are two so-called covariants, which we denote g_{4} and g_{6}:

(3.6.4)  \begin{aligned}g_{4}(X, Y) & = (3b^{2} − 8ac)X^{4} + 4(bc − 6ad)X^{3}Y + 2(2c^{2} − 24ae − 3bd)X^{2}Y^{2} \\ & + 4(cd − 6be)XY^{3} + (3d^{2} − 8ce)Y^{4}, \\ g_{6}(X, Y) & = (b^{3} + 8a^{2}d − 4abc)X^{6} + 2(16a^{2}e + 2abd − 4ac^{2} + b^{2}c)X^{5}Y \\ & + 5(8abe + b^{2}d − 4acd)X^{4}Y^{2} + 20(b^{2}e − ad^{2})X^{3}Y^{3} \\ & − 5(8ade + bd^{2} − 4bce)X^{2}Y^{4} − 2(16ae^{2} + 2bde − 4c^{2}e + cd^{2})XY^{5} \\ & − (d^{3} + 8be^{2} − 4cde)Y^{6}. \end{aligned}

Wherever convenient, we will also denote by g the homogenized polynomial g(X, Y) = aX^{4} + bX^{3}Y + cX^{2}Y^{2} + dXY^{3} + eY^{4}. These three homogeneous polynomials satisfy an algebraic identity, or syzygy:

(3.6.5)  27g_{6}^{2} = g_{4}^{3} − 48Ig^{2}g_{4} − 64Jg^{3}.

Later we will also need a simpler form of this syzygy; set

(3.6.6)  p = g_{4}(1, 0) = 3b^{2} − 8ac  and  r = g_{6}(1, 0) = b^{3} + 8a^{2}d − 4abc;

these quantities are called seminvariants of g. Substituting (X, Y) = (1, 0) in the covariant syzygy (3.6.5) gives an identity (the seminvariant syzygy) between these seminvariants:

(3.6.7)  27r^{2} = p^{3} − 48Ia^{2}p − 64Ja^{3}.

We will make use of this equation in our search for quartics with given invariants, where it will allow us to set up a quadratic sieve.

It follows from the covariant syzygy (3.6.5), by simple substitution, that the map

(3.6.8)  ξ : (x, y) ↦ (\frac{3g_{4}(x, 1)}{(2y)^{2}}, \frac{27g_{6}(x, 1)}{(2y)^{3}})

maps rational points (x, y) on H (satisfying y^{2} = g(x, 1)) to rational points on E_{I,J}, thus defining a rational map ξ, of degree 4, from H(ℚ) to E_{I,J}(ℚ). We are using affine coordinates

here; the points at infinity on H map to (\frac{3p}{4a}, \frac{±27r}{(4a)^{3/2}} ), which are rational if and only if a is a square.

We now have the following facts (see [20] for details):

• If R ∈ H(ℚ) with P = ξ(R) ∈ E_{I,J}(ℚ), then the coset of P modulo 2E_{I,J}(ℚ) is independent of R, and of the particular quartic g up to equivalence; in fact, equivalences between quartics induce rational maps between the associated homogeneous spaces, and the covariant property of g_{4} and g_{6} ensures that corresponding rational points on the homogeneous spaces have the same image in E_{I,J}(ℚ).

• Each rational point P = (x, y) ∈ E_{I,J}(ℚ) arises as the image of a rational point on some quartic g with invariants I and J: explicitly, one can take the rational point at infinity on the quartic with coefficients (a, b, c, d, e) = (1, 0, −x/6, y/27, I/12 − x^{2}/432); the equivalence class of g depends only on the coset of P modulo 2E_{I,J}(ℚ).

• The equivalence classes of everywhere locally soluble quartics with invariants I and J form a finite elementary abelian 2-group, isomorphic to the 2-Selmer group S^{(2)}(E/ℚ).

• The equivalence classes of soluble quartics with invariants I and J form a finite elementary abelian 2-group isomorphic to E(ℚ)/2E(ℚ); the identity is the trivial class, consisting of quartics with a rational root.

• More generally, when E has no 2-torsion, for any extension field K of ℚ there is a bijection between the roots of g(x) in K and the solutions Q ∈ E_{I,J}(K) to the equation 2Q = P (where P = ξ(R) for R ∈ H(ℚ) as above). In particular, non-trivial quartics are irreducible in this case. We will use this fact with K = ℝ later.

We therefore classify the set of equivalence classes of quartics with invariants I and J as follows:

(0) the trivial class consists of those quartics g(x) which have a rational root. These are elliptic curves isomorphic to E over ℚ.

(1) those which have a rational point: these are also elliptic curves, isomorphic to E over ℚ.

(2) those which have points everywhere locally.

(3) those which fail to have points everywhere locally.

Let the number of inequivalent quartics in the first three sets be n_{0} = 1, n_{1} and n_{2}. (Those in the last set will not be used.) Because of the group structure, each of these numbers is a power of 2. We write n_{i} = 2^{e_{i}} for i = 1, 2.

As in the case of descent via 2-isogeny, Galois cohomology gives an exact sequence

0 → E(ℚ)/2E(ℚ) → S^{(2)}(E/ℚ) → � (E/ℚ)[2] → 0.

Thus the quotient of S^{(2)}(E/ℚ) by the image of E(ℚ) is isomorphic to � (E/ℚ)[2], the 2torsion subgroup of the Tate–Shafarevich group � (E/ℚ). So it is the points of order 2 in � (E/ℚ), if any, which account for the possible existence of homogeneous spaces which have points everywhere locally but not globally, and we have

|� (E/ℚ)[2]| = n_{2}/n_{1}.

As before, the potential practical difficulty lies in determining whether each homogeneous space H has a rational point, as there is no known algorithm to do this in general. Again, for the vast majority of the curves in the tables, we found a rational point easily on each space which was everywhere locally soluble, which not only determined the rank of E, but also implied that the Tate–Shafarevich group had no 2-torsion. The only example with n_{1} < n_{2} in the tables (for a curve with no 2-torsion) is curve 571A1, where n_{1} = 1 and n_{2} = 4; here the rank is 0, and |� (E/ℚ)[2]| = 4; the Birch–Swinnerton-Dyer conjecture predicts |� (E/ℚ)| = 4.

The steps of the algorithm are as follows: first we determine the pair or pairs of integral invariants (I, J) such that every quartic associated with our curve E is equivalent to one with integer coefficients and these invariants. There will be either one or two such pairs. For each pair (I, J), we find a finite set of quartics with invariants (I, J) such that every non-trivial, everywhere locally soluble quartic with these invariants is equivalent to one in the list. This is the most time-consuming step, as the search region can be very large when I and J are large. Now we must test the quartics in our list pairwise for equivalence, discarding those equivalent to earlier ones; look for rational points; and test everywhere local solubility. Again, there may be quartics where we do not find rational points despite their having points everywhere locally, so that although we can always (given enough time) determine n_{2}, we may in some cases only find bounds on n_{1}. Since n_{1} = |E(ℚ)/2E(ℚ)|, we can then compute the rank r, or bounds on the rank. Usually, E will have no rational 2-torsion, or we would probably be using descent via 2-isogeny, and then simply 2^{r} = n_{1}.

We now consider each of these steps in more detail.

Step 1: Determining the invariants (I, J).

Given an integral quartic g with invariants I and J, we must consider the question of whether there exists an equivalent integral quartic with smaller invariants. The smaller invariants will have the form λ^{−4}I, λ^{−6}J with λ ∈ ℚ^{∗}. In [3, Lemmas 3–5], conditions are stated under which g is equivalent to an integral quartic with invariants p^{−4}I, p^{−6}J for a prime p; we call such a quartic p-reducible, otherwise p-minimal. Clearly a necessary condition for reducibility is that p^{4} | I and p^{6} | J. We say that the pair (I, J) is p-reducible if every integral quartic with these invariants which is p-adically soluble is equivalent to an integral quartic with invariants p^{−4}I and p^{−6}J.

The question of p-reducibility is almost completely settled by the following proposition. The result is simplest for primes greater than 3, but even for these it is important to note that the assumption of p-adic solubility is necessary for reduction to be possible when the divisibility conditions are satisfied.

Proposition 3.6.1. Let I and J be integers such that Δ = 4I^{3} − J^{2} ≠ 0.

(1) If p is a prime and p ≥ 5, then (I, J) is p-reducible if and only if p^{4} | I and p^{6} | J.

(2) (I, J) is 3-reducible if and only if either 3^{5} | I and 3^{9} | J, or 3^{4} || I, 3^{6} || J and 3^{15} | Δ.

(3) (I, J) is 2-reducible if 2^{6} | I, 2^{9} | J and 2^{10} | 8I + J.

This proposition is stated in [3] as Lemmas 3–5, but only the proof of Lemma 3 (covering the case p ≥ 5) is given there. Complete proofs in all cases (which are elementary though somewhat lengthy) can be found in [52].

Note that for p = 2 we only have sufficient conditions for reducibility. Because of this, we will sometimes have to consider two pairs of invariants, a smaller pair (I_{0}, J_{0}) and a larger pair (16I_{0}, 64J_{0}). However, when searching for integral quartics with the larger invariants, we may assume that the quartic cannot be 2-reduced, and this provides us with useful congruence conditions on the coefficients of such a quartic. We state these here.

Proposition 3.6.2. Let g be an integral 2-adically soluble quartic whose invariants satisfy 2^{4} | I and 2^{6} | J, such that

(1) g is not equivalent to an integral quartic with invariants 2^{−4}I and 2^{−6}J;

(2) g is not equivalent to an integral quartic with the same invariants I and J and smaller

leading coefficient a. Then the coefficients of g satisfy

(a) 2 ∤ a, 2^{2} | b, 2 | c, 2^{4} ∤ e and 2^{4} ∤ a + b + c + d + e; or

(b) 2 || a, 2^{2} | b, 2^{2} | c, 2^{3} ∤ e and 2^{3} ∤ a + b + c + d + e. Moreover, if 2^{6} | I and 2^{7} | J, then we must have

(a’) 2 ∤ a, 2^{2} | b, 2^{2} || c, 2^{3} | d, and 2^{2} || e; or (b’) 2 ∤ a, 2^{2} | b, 2^{2} || c − 2a + 3b, 2^{3} | d − b and 2^{2} || a + c + e.

The first set of conditions stated here were given in [3]; the second set are from [52], which contains complete proofs in both cases.

Using this proposition, we may ensure that few of the quartics we find when searching the larger pair of invariants are equivalent to one with smaller invariants. More significantly in terms of running time, we have extra congruence conditions to apply when searching for the larger invariants, which speeds up this search.

It would appear that rational points in E(ℚ) whose quartics have the larger pair of invariants lie in certain components of the 2-adic locus E(ℚ_{2}). Further study of this would be very useful, since if the search for quartics with the larger pair of invariants could be eliminated or curtailed, it could result in a major saving of time in the algorithm.

In practice, suppose that our original curve E is given by a minimal equation, with invariants c_{4} and c_{6}. We set I = c_{4} and J = 2c_{6}. Clearly the pair (I, J) is p-minimal for p ≥ 5: for if p^{4} | I and p^{6} | J then p^{−4}c_{4} and p^{−6}c_{6} would be integral invariants of an elliptic curve, contradicting minimality of E, and similarly the pair (p^{4}I, p^{6}J) is certainly p-reducible by Proposition 3.6.1(1). Less obvious is that (I, J) is also 3-minimal; using Kraus’s conditions, it is easy to check first that (3^{4}I, 3^{6}J) is certainly 3-reducible (one needs here that ord_{3}(c_{6}) ≠ 2), and then that (I, J) itself is not 3-reducible, using Proposition 3.6.1(2).

For p = 2, the best we can do is the following. First set I = c_{4} and J = 2c_{6}. Replace (I, J) by (2^{−4}I, 2^{−6}J) if 2^{4} | I and 2^{6} | J; the resulting pair (I, J) (which will not be further divisible by 2) will be the basic pair of invariants. Then we also use the pair (16I, 64J) unless 4 | I, 8 | J and 16 | (2I + J).

The result of this step is then to produce either one or two pairs of invariants (I, J). In the latter case, the following steps must be carried out with both pairs separately.

Step 2: Finding the quartics with given I and J.

We now have a fixed pair of invariants (I, J) with Δ = 4I^{3} − J^{2} ≠ 0, and we wish to find all integral quartics with these invariants, up to equivalence. We classify the quartics g(x) into types, according as g(x) has no real roots (type 1), four real roots (type 2) or two real roots (type 3). When Δ < 0 only type 3 is possible, while if Δ > 0, only types 1 and 2 are possible. For each relevant type, we now determine a finite list of quartics of that type with the given invariants such that every soluble quartic with these invariants is equivalent to at least one on the list. We can ignore quartics which are negative definite (type 1 with a < 0), since they will not be soluble over ℝ. For each type, we will determine a finite region of (a, b, c)-space such that every quartic with invariants I and J is equivalent to at least one in this region.

As observed above, the number of real roots of g(x) is equal to the number of points Q ∈ E(ℝ) satisfying 2Q = P, where P ∈ E(ℝ) is the image under the map ξ of any real point on the homogeneous space H with equation y^{2} = g(x). When Δ < 0, the real locus is in one component, and E(ℝ) is isomorphic to the circle group, which is 2-divisible with two 2-torsion points, so in this case the equation 2Q = P has exactly two solutions for all P ∈ E(ℝ). This agrees with the observation just made, that quartics with negative discriminant Δ will all have exactly two real roots.

Consider further the case Δ > 0. Now E(ℝ) has two components, the connected component of the identity E^{0}(ℝ) and a second component which we call the ‘egg’. There are four 2-torsion points, and 2E(ℝ) = E^{0}(ℝ). There are therefore two possibilities for a point P ∈ E(ℝ) and its associated real quartic: if P ∈ E^{0}(ℝ), then there are four solutions Q to 2Q = P, and P will be associated to a quartic of type 2 with four real roots. On the other hand, if P ∉ E^{0}(ℝ), then there are no solutions and the quartic associated to P will be of type 1, with no real roots.

The image of E(ℚ) in E(ℝ)/2E(ℝ) has order 2 or 1, depending on whether or not there are any rational points on the egg. Thus there are two sub-cases to the case Δ > 0: if E(ℚ) ⊂ E^{0}(ℝ), then there are no rational points on the egg, the index is 1, and there will be no soluble quartics of type 1; on the other hand, if E(ℚ) ⊄ E^{0}(ℝ), then there are rational points on the egg, the index is 2, and there are equal numbers of (equivalence classes of) soluble quartics of types 1 and 2. Those of type 2 will lead to rational points on E(ℚ) ∩ E^{0}(ℝ), while those of type 1 will lead to rational points on the egg.

To take advantage of this in practice, when Δ > 0 we will first look for quartics of type 2; let the number of these be n_{1}^{+}, where n_{1}/n_{1}^{+} is either 1 or 2. At this stage we will already know the rank to within one, since if we set r^{+} = log_{2}(n_{1}^{+}) then (assuming no rational 2-torsion) we have either r = r^{+} or r = r^{+}+1. Then we start to look for quartics of type 1; as soon as we find one which is soluble, then we may abort the search for type 1 quartics at that point, and assert that r = r^{+} + 1. On the other hand, if we complete the search for quartics of type 1 without finding any soluble ones, then we will know that r = r^{+}, and we will have proved that there are no rational points on the egg. An example of the second possibility happens with the curve E = [0, 0, 1, −529, −3042] (which is the −23-twist of the curve [0, 0, 1, −1, 0] with conductor 37 and rank 1), which has rank 1 with generator (46, 264) on the identity component, and no rational points on the egg.^{8}

If we happened to know in advance that there were rational points on the egg (perhaps by a short preliminary search for such points with small height), then we would already know that r = r^{+} + 1, and we would not need to search for type 1 quartics at all.

In order to find all integral quartics of a given type (up to equivalence) we proceed as follows. First, following [3], we determine bounds on the coefficients a, b and c. We also set up a sieve based on the seminvariant syzygy (3.6.7) to speed up our search through this region of (a, b, c)-space. For triples (a, b, c) in the region which pass the sieve, we solve for d and e and ensure that they are integral. Finally, we check that the quartic we have constructed satisfies any further congruence conditions we require (for example, when we are using the larger pair of invariants).

8Thanks to Nelson Stephens for this example.

The method for bounding the coefficients which is developed in [3] involves using the auxiliary (resolvent) cubic equation

(3.6.9)  ϕ^{3} − 3Iϕ + J = 0

which will have one real root (type 3) or three real roots (types 1 and 2), since its discriminant is 27Δ. Indeed, ϕ is a root of (3.6.9) if and only if (−3ϕ, 0) is a point of order 2 on the curve E_{I,J}.

In each case, the bound for b arises simply from the fact that the quartics g(x) and g(x + k) are equivalent, and the coefficients of the latter are (a, b + 4ak, …), so that we may assume that b is reduced modulo 4a. Also, note that the bounds on c are effectively bounds on the seminvariant 8ac − 3b^{2} = −p, which is how they arise in [3].

Bounds for (a, b, c): Type 1. Here we may assume a > 0 for real solubility. Order the three real roots of (3.6.9) as ϕ_{1} > ϕ_{2} > ϕ_{3}, and set K = (4I − ϕ_{1}^{2})/3. Then the bounds on a, b, c are

0 < a ≤ \frac{K + K^{\frac{1}{2}} ϕ_{1}}{3K^{\frac{1}{2}} + ϕ_{1} + 2ϕ_{2}} ;
−2a < b ≤ 2a;
\frac{4aϕ_{2} + 3b^{2}}{8a} ≤ c ≤ \frac{4aϕ_{1} + 3b^{2}}{8a} .

Bounds for (a, b, c): Type 2. This subdivides into subtypes according as a > 0 or a < 0. For a > 0 we take ϕ_{1} > ϕ_{2} > ϕ_{3} and search the region

0 < a ≤ \frac{I − ϕ_{2}^{2}}{3(ϕ_{2} − ϕ_{3})};
−2a < b ≤ 2a;
\frac{4aϕ_{2} − \frac{4}{3}(I − ϕ_{2}^{2}) + 3b^{2}}{8a} ≤ c ≤ \frac{4aϕ_{3} + 3b^{2}}{8a} .

Then for a < 0 we take ϕ_{1} < ϕ_{2} < ϕ_{3} and search over

0 < − a ≤ \frac{I − ϕ_{2}^{2}}{3(ϕ_{3} − ϕ_{2})};
−2|a| < b ≤ 2|a|;
\frac{4aϕ_{2} − \frac{4}{3}(I − ϕ_{2}^{2}) + 3b^{2}}{8a} ≥ c ≥ \frac{4aϕ_{3} + 3b^{2}}{8a} .

Bounds for (a, b, c): Type 3. Here we let ϕ be the unique real root of (3.6.9), and search

\frac{1}{3}ϕ − \sqrt{\frac{4}{27}(ϕ^{2} − I)} ≤ a ≤ \frac{1}{3}ϕ + \sqrt{\frac{4}{27}(ϕ^{2} − I)};
−2|a| < b ≤ 2|a|;
\frac{9a^{2} − 2aϕ + \frac{1}{3}(4I − ϕ^{2}) + 3b^{2}}{8|a|} ≤ c.sign(a) ≤ \frac{4aϕ + 3b^{2}}{8|a|} .

The syzygy sieve. Recall the seminvariant syzygy

(3.6.7)  27r^{2} = p^{3} − 48Ia^{2}p − 64Ja^{3} = s(a, p),

say, where p = 3b^{2} − 8ac and r = b^{3} + 8a^{2}d − 4abc. For fixed I, J the expression s(a, p) is a polynomial in a, b and c, which we require to be 27 times an integer square. We can set up a quadratic sieve as follows: for each of several sieving moduli m we create and initialize an m × m array indicating whether s(a, p) is 27 times a square modulo m, for each pair (a, p) modulo m. We take one of the moduli to be 9 and use it to force the right-hand side of (3.6.7) to be divisible by 27; it will certainly be positive, as this is ensured by the bounds on c.

For each (a, b, c) in the region searched, we check that it passes the sieving test; it is then quite likely that s(a, p) will be 27 times a square, since it is so modulo a large modulus and is positive. We then test whether this is the case, discarding (a, b, c) if not, and if so we then find r. We can take r > 0, since the quartics with coefficients (a, b, c, d, e) and (a, −b, c, −d, e) are equivalent, with opposite signs of their respective r-seminvariants. In fact, we treat the triples (a, ±b, c) together in practice.

Implementation note: It is worth pointing out that a large proportion of the running time of our algorithm is spent testing whether large integers are squares (given that they are positive and congruent to squares modulo several carefully chosen moduli), and find their integer square root if so. This is needed here, and in our searches for rational points, both on the elliptic curve directly, and on the homogeneous spaces. Hence it is crucial that we have access to efficient procedures for this in the multiprecision integer package we use.

Solving for d and e. Given integers a, b, c, r satisfying (3.6.7) with p = 3b^{2} − 8ac, we can solve for d and e, setting

d = (r − b^{3} + 4abc)/(8a^{2})  and  e = (I + 3bd − c^{2})/(12a).

This will certainly give rational values for d and e; we must check that they are integral, discarding the triple (a, b, c) if not. If they are, we have integral coefficients (a, b, c, d, e) of a quartic g(x) with invariants I and J in the search region, which we add to our list for further processing.

Solving for the roots of g(x). For later use, when we check for triviality, and again when we search for rational points on the homogeneous spaces, we will need to know the real roots of the quartic g(x) we have constructed. Although the formulae for finding the roots of quartic are well-known, we give them here: since we already know the roots of the resolvent cubic, there is very little work remaining.

For i = 1, 2, 3 we set z_{i} = (4aϕ_{i} + p)/3 where the ϕ_{i} are the three roots of (3.6.9). The product of these quantities is r^{2} (from (3.6.7) again), and we form their square roots with product r by setting w_{1} = \sqrt{z_{1}}, w_{2} = \sqrt{z_{2}}, and w_{3} = r/(w_{1}w_{2}). Then the roots of g(x) are

x_{1} = ( w_{1} + w_{2} − w_{3} − b)/(4a),
x_{2} = ( w_{1} − w_{2} + w_{3} − b)/(4a),
x_{3} = (−w_{1} + w_{2} + w_{3} − b)/(4a),
x_{4} = (−w_{1} − w_{2} − w_{3} − b)/(4a).

We will not give here a pseudo-code algorithm for the search for quartics, as it is straightforward in principle, although in practice it needs careful book-keeping. As this is the most time-consuming part of the whole procedure, particularly when the second, larger, pair of invariants must be used, it is important to make the implementation code as efficient as possible.

At the end of this step we will have a list of quartics with the desired invariants. We now discard any which are equivalent to earlier ones, or are not locally soluble at some prime p, and try to find rational roots on the remainder. In practice we may choose to apply these tests in a different order, such as not bothering to check equivalences between quartics which are not locally soluble.

Step 3: Testing triviality.

For each quartic g(x) in the list, we already know its roots x to reasonable precision. If x is rational, then ax is integral, which we can test. If we suspect that ax is equal to an integer n to within some working tolerance, we can check whether n/a is a root of g(x) using exact arithmetic.

Step 3: Testing equivalence of quartics.

With each quartic we find with the right invariants, we store its coefficients, type, roots and seminvariants p and r. We also compute and store the number of roots of the quartic (including roots at infinity) modulo each of several primes not dividing its discriminant, as these numbers are clearly invariant under equivalence.^{9}

When testing equivalence of two quartics, we first check that their invariants and type are the same, as well as their numbers of roots modulo these primes. If this is the case, we use a general test for equivalence (valid over any field) from [20], which we state here.^{10}

Proposition 3.6.3. Let g_{1} and g_{2} be quartics over the field K, both having the same invariants I and J, and with leading coefficients a_{i} and seminvariants p_{i} and r_{i} for i = 1, 2. Then g_{1} is equivalent to g_{2} over K if and only if the quartic u^{4} − 2pu^{2} − 8ru + s has a root in K, where

p = (32a_{1}a_{2}I + p_{1}p_{2})/3,
r = r_{1}r_{2},  and
s = (64I(a_{1}^{2}p_{2}^{2} + a_{2}^{2}p_{1}^{2} + a_{1}a_{2}p_{1}p_{2}) − 256a_{1}a_{2}J(a_{1}p_{2} + a_{2}p_{1}) − p_{1}^{2}p_{2}^{2})/27.

The quantities p, r and s in this proposition will be integers when g_{1} and g_{2} are integral. Converting the proposition into an algorithm is straightforward.

Step 5: Testing local and global solubility.

This is carried out using the procedures and strategy described earlier.

Step 6: Final computation of the rank.

The number of quartics found (up to equivalence) which are everywhere locally soluble is n_{2}, the order of the 2-Selmer group. This must be a power of 2, say n_{2} = 2^{e_{2}}, which serves as a check on our procedures. The number n_{1} with a rational point is also a power of 2, say n_{1} = 2^{e_{1}}, equal to the order of E(ℚ)/2E(ℚ). If we have found rational points on all n_{2} locally soluble quartics, then certainly n_{1} = n_{2}, so that � (E/ℚ)[2] is trivial, and the rank of E(ℚ) is e_{1} − e_{0} where |E(ℚ)[2]| = 2^{e_{0}} with e_{0} = 0, 1 or 2. The rank is equal to the Selmer rank e_{2} − e_{0} in this case. (Usually e_{0} = 0 when we are using this method.)

As before, we may not have found global points on all the locally soluble quartics; if the number on which we have points is ñ_{1} with ñ_{1} < n_{2} then we only know that ñ_{1} ≤ n_{1} ≤ n_{2}. If ñ_{1} is not a power of 2, we will know that n_{1} > ñ_{1}, so that at least some of our locally soluble quartics must have rational points which we have not found. In this case, we replace ñ_{1} by the next highest power of 2, say ñ_{1} = 2ẽ_{1}. Then we have bounds on the rank, namely

9This was suggested to us by S. Siksek.

10The algorithm presented here only applies to quartics. In the First Edition we presented a different algorithm, described in [3], which is messier to implement, but which generalizes more readily to more general situations, such as testing the equivalence of binary forms of higher degree.

ẽ_{1} − e_{0} ≤ e_{1} − e_{0} = rank(E(ℚ)) ≤ e_{2} − e_{0},

and on the order of � (E/ℚ)[2]:

|� (E/ℚ)[2]| ≤ n_{2}/ñ_{1}.

One final point: from the Selmer conjecture, we expect the Selmer rank e_{2} − e_{0} to differ from the actual rank e_{1} − e_{0} by an even number, so that e_{2} ≡ e_{1} (mod 2). This would also follow from the conjecture that � (E/ℚ) is finite, since then its order is known to be a perfect square, so that n_{2}/n_{1} must be a square. So if we find that e_{2} ≢ ẽ_{1} (mod 2), then we suspect that the rank is at least one more than our lower bound, and can output a comment to this effect, though of course we will not have proved that the rank is greater than our lower bound. In some cases, such as for a modular curve where we know the sign of the functional equation, we may have other conjectural evidence for the parity of the rank.

Step 7: Recovering points on E.

Each quartic g(x) for which the homogeneous space y^{2} = g(x) has a rational point R leads to a rational point P = ξ(R) on the model E_{I,J} of our curve E, via the formula (3.6.8) given above. If we apply this formula to all the inequivalent quartics with rational points which we found in computing the rank of E, we will have a complete set of coset representatives for 2E(ℚ) in E(ℚ), provided that ñ_{1} = n_{1}. In cases where we have rounded up ñ_{1} to the nearest power of 2, we will still have generators for E(ℚ)/2E(ℚ), and can fill in the missing coset representatives if we wish.

This completes our description of algorithms for determining the Mordell-Weil group E(ℚ).

3.7 The period lattice

In this section we show how to compute the complex periods for an elliptic curve defined over the complex numbers. We used this in our investigation of modular curves to check that the exact integral equations we found (after rounding the approximate computed values of c_{4} and c_{6}) did have the correct periods; and also in our method for computing isogenous curves, which we describe in the following section.

Let E be an elliptic curve defined over the complex numbers ℂ, given by a Weierstrass equation. We wish to compute periods λ_{1} and λ_{2} which are a ℤ-basis for the period lattice Λ of E. We do this using Gauss’s arithmetic–geometric mean (AGM) algorithm. Write the equation for E in the form

(y + \frac{a_{1}x + a_{3}}{2} )^{2} = x^{3} + \frac{b_{2}}{4} x^{2} + \frac{b_{4}}{2} x + \frac{b_{6}}{4} = (x − e_{1})(x − e_{2})(x − e_{3}),

where the roots e_{i} are found as complex floating-point approximations (using Cardano’s formula, say). Then the periods are given by

(3.7.1)  \begin{aligned}λ_{1} & = \frac{π}{AGM(\sqrt{e_{3} − e_{1}}, \sqrt{e_{3} − e_{2}})}, \\ λ_{2} & = \frac{πi}{AGM(\sqrt{e_{3} − e_{1}}, \sqrt{e_{2} − e_{1}})}. \end{aligned}

Notice that in general this involves the AGM of pairs of complex numbers. This is a multivalued function: at each stage of the AGM algorithm we replace the pair (z, w) by (\sqrt{zw}, \frac{1}{2}(z + w)), and must make a choice of complex square root. It follows from work of Cox (see [11]) that while a different set of choices does lead to a different value for the AGM, the periods we obtain this way will nevertheless always be a ℤ-basis for the full period lattice Λ. We have found this to be the case in practice, where we always choose a square root with positive real part, or with positive imaginary part when the real part is zero. The computation of λ_{1} and λ_{2} by this method is very fast, as the AGM algorithm converges extremely quickly, even in its complex form. As a check on the values obtained, in each case we recomputed the invariants c_{4} and c_{6} of each curve from these computed periods λ_{1} and λ_{2}, using the standard formulae given in Chapter II; in every case we obtained the correct values (known exactly from the coefficients of the minimal Weierstrass equation) to within computational accuracy.

If the curve is defined over ℝ, we can avoid the use of the complex AGM, and also arrange that λ_{1} is a positive real period, as follows. First suppose that all three roots e_{i} are real; order the roots so that e_{3} > e_{2} > e_{1}, and take the positive square root in the above formulae. Then we may use the usual AGM of positive reals in (3.7.1), and thus obtain a positive real value for λ_{1} and a pure imaginary value for λ_{2}. This is the case where the discriminant Δ > 0 and the period lattice is rectangular. When Δ < 0 there is one real root, say e_{3}, and e_{2} = \overline{e_{1}}. If \sqrt{e_{3} − e_{1}} = z = s + it with s > 0 then \sqrt{e_{3} − e_{2}} = \overline{z} = s − it, so that λ_{1} = π/AGM(z, \overline{z}) = π/AGM(|z|, s) which is also real and positive.

3.8 Finding isogenous curves

Given an elliptic curve E defined over ℚ, we now wish to find all curves E′ isogenous to E over ℚ. The set of all such curves is finite (up to isomorphism), and any two curves in the isogeny class are linked by a chain of isogenies of prime degree l. Thus it suffices to be able to compute l-isogenies for prime l, if we can determine those l for which rational l-isogenies exist. The latter question can be rather delicate in general, and we have to have a completely automatic algorithmic procedure if we are to apply it to several thousand curves, such as we had to when preparing the tables.

When the conductor N of E is square-free, so that E has good or multiplicative reduction at all primes, E is called semi-stable. In this case, a result of Serre (see [53]) says that either E or the isogenous curve E′ has a rational point of order l, and so by Mazur’s result already mentioned, l can only be 2, 3, 5 or 7. Moreover, if a curve E possesses a rational point of order l, then the congruence 1 + p − a_{p} ≡ 0 (mod l) holds for all primes p not dividing Nl, so the presence of such a point is easy to determine, even if it is not E itself but the isogenous curve E′ which possesses the rational l-torsion, since the trace of Frobenius a_{p} is isogeny-invariant.

If E is not semi-stable we argue as follows. The existence of a rational l-isogeny is purely a function of the j-invariant j of E: in fact, pairs (E, E′) of l-isogenous curves parametrize the modular curve X_{0}(l) whose non-cuspidal points are given by the pairs (j(E), j(E′)). For l = 2, 3, 5, 7 or 13 the genus of X_{0}(l) is zero, and infinitely many rational j occur. The only other values of l for which rational l-isogenies occur are l = 11, 17, 19, 37, 43, 67, and 163, and these occur for only a small finite number of j-invariants (see below). The fact that no other l occur is a theorem of Mazur (see [39] and [40]), related to the theorem limiting the rational torsion which we quoted earlier in Section 3.3 of this chapter. These extra values occur only for curves with CM (complex multiplication, see the next section), apart from l = 17 (where X_{0}(l) has genus 1) and the exotic case l = 37 studied by Mazur and Swinnerton–Dyer in [41] (where X_{0}(l) has genus 2).

For isogenies of non-prime degree m, the degrees which occur are: m ≤ 10, and m = 12, 16 18, and 25 (where X_{0}(l) has genus 0, infinitely many cases); and finally m = 14, 15, 21, and 27. The latter occur first for conductors N = 49 (with CM), N = 50, N = 162 and N = 27 (with CM) respectively. See [2, pages 78–80] for more details.

Thus our procedure is:

• If N is square-free, try l = 2, 3, 5, 7 only;

• else try l = 2, 3, 5, 7 and 13 in all cases; and

• if j(E) = −2^{15}, −11^{2}, or −11 · 131^{3}, try also l = 11;

• if j(E) = −17^{2} · 101^{3}/2 or −17 · 373^{3}/2^{17}, try also l = 17;

• if j(E) = −96^{3}, try also l = 19;

• if j(E) = −7 · 11^{3} or −7 · 137^{3} · 2083^{3}, try also l = 37;

• if j(E) = −960^{3}, try also l = 43;

• if j(E) = −5280^{3}, try also l = 67;

• if j(E) = −640320^{3}, try also l = 163.

Now we turn to the question of finding all curves (if any) which are l-isogenous to our given curve E for a specific prime l. The kernel of the isogeny is a subgroup A of E(\overline{ℚ}) which is defined over ℚ, but the points of A may not be individually rational points. If we have the coordinates of the points of a subgroup of E of order l defined over K, we may use Vélu’s formulae in [68] to find the corresponding l-isogenous curve. Finding such coordinates by algebraic means is troublesome, except when the subgroup is point-wise defined over K, and instead we resort to a floating-point method.

The case l = 2 is simpler to describe separately. Obviously in this case the subgroup of order 2 defined over ℚ must consist of a single rational point P of order 2 together with the identity. We have already found such points, if any, in computing the torsion. There will be 0, 1 or 3 of them according to the number of rational roots of the cubic 4x^{3} + b_{2}x^{2} + 2b_{4}x + b_{6}. If x_{1} is such a root, then P = (x_{1}, y_{1}) has order 2, where y_{1} = −(a_{1}x_{1} + a_{3})/2. As a special case of Vélu’s formulae we find that the isogenous curve E′ has coefficients [a′_{1}, a′_{2}, a′_{3}, a′_{4}, a′_{6}] = [a_{1}, a_{2}, a_{3}, a_{4} − 5t, a_{6} − b_{2}t − 7w] where

t = (6x_{1}^{2} + b_{2}x_{1} + b_{4})/2  and  w = x_{1}t.

Note that the point (x_{1}, y_{1}) need not be integral even when E has integral coefficients a_{i}, but that 4x_{1} and 8y_{1} are certainly integral, by the formulae given; thus the model just given for the isogenous curve may need scaling by a factor of 2 to make it integral.

The simpler formula for a curve in the form y^{2} = x^{3} + cx^{2} + dx and the point P = (0, 0) was given in the previous section: the formulae just given take the curve [0, c, 0, d, 0] to [0, c, 0, −4d, −4cd], which transforms to [0, −2c, 0, c^{2} − 4d, 0] after replacing x by x − c. The relation between the two formulae is given by c = 12x_{1} + b_{2} and d = 16t.

For reference we give here similar algebraic formulae for l-isogenies for l = 3 and l = 5, from Laska’s book [35]. In each case we assume that the curve E is given by an equation of the form y^{2} = x^{3} + ax + b, and the isogenous curve E′ by y^{2} = x^{3} + Ax + B. Each subgroup of E of order l is determined by a rational factor of degree (l − 1)/2 of the l-division polynomial of degree (l^{2} − 1)/2, whose roots are the x-coordinates of the points in the subgroup. The simplest case is l = 3, where there is just one x-coordinate, which must be rational.

l = 3. Let ξ be a root of the 3-division polynomial 3x^{4} + 6ax^{2} + 12bx − a^{2}. Then the 3-isogenous curve E′ is given by

A = −3(3a + 10ξ^{2})
B = −(70ξ^{3} + 42aξ + 27b).

l = 5. Let x^{2} + h_{1}x + h_{2} be a rational factor of the 5-division polynomial 5x^{12} + 62ax^{10} + 380bx^{9} − 105a^{2}x^{8} + 240abx^{7} − (300a^{3} + 240b^{2})x^{6} − 696a^{2}bx^{5} − (125a^{4} + 1920ab^{2})x^{4} − (1600b^{3} + 80a^{3}b)x^{3}−(50a^{5}+240a^{2}b^{2})x^{2}−(100a^{4}b+640ab^{3})x+(a^{6}−32a^{3}b^{2}−256b^{4}). Then the 5-isogenous curve E′ is given by

A = −19a − 30(h_{1}^{2} − 2h_{2})
B = −55b − 14(15h_{1}h_{2} − 5h_{1}^{3} − 3ah_{1}).

A similar formula is given in [35] for l = 7, where A and B are given in terms of a, b and the coefficients of a factor x^{3} + h_{1}x^{2} + h_{2}x + h_{3} of the 7-division polynomial. Rather than take up space by giving the latter here, we refer the reader to [35, page 72].

Now we turn to Vélu’s formulae in the case of an odd prime l. Let P = (x_{1}, y_{1}) be a point of order l in E(\overline{ℚ}), and set kP = (x_{k}, y_{k}) for 1 ≤ k ≤ (l − 1)/2. Define

t_{k} = 6x_{k}^{2} + b_{2}x_{k} + b_{4}  and  u_{k} = 4x_{k}^{3} + b_{2}x_{k}^{2} + 2b_{4}x_{k} + b_{6},

and then set

t = ∑_{k=1}^{(l−1)/2} t_{k}  and  w = ∑_{k=1}^{(l−1)/2} (u_{k} + x_{k}t_{k}) .

Then the isogenous curve E′ has coefficients [a_{1}, a_{2}, a_{3}, a_{4} − 5t, a_{6} − b_{2}t − 7w] as before. Again, these may not be integral, even when the original coefficients were; but since the x_{k} are the roots of a polynomial of degree (l − 1)/2 with integral coefficients and leading coefficient l^{2} (the so-called l-division equation), we must have l^{2}x_{k} integral. Thus a scaling factor of l will certainly produce an integral equation.

We make these remarks on integrality as our method is to find the coordinates x_{k} and y_{k} as real floating-point approximations, and thus to determine the coefficients of any curves l-isogenous to E over ℝ; there will always be exactly two such curves over ℝ, but of course they will not necessarily be defined over ℚ. As we will only know the coefficients a′_{i} of the isogenous curves approximately, we wish to ensure that if they are rational then they will in fact be integral, so that we will be able to recognize them as such.

First we find the period lattice Λ of E, as described in the previous section. The ℤ-basis [λ_{1}, λ_{2}] of Λ is normalized as follows: there are two cases to consider, according as Δ > 0 (first or ‘harmonic’ case) or Δ < 0 (second or ‘anharmonic’ case). In both cases λ_{1} is real (the least positive real period); in the first case, λ_{2} is pure imaginary, while in the second case, 2λ_{1} − λ_{2} is pure imaginary. We can also ensure that τ = λ_{2}/λ_{1} is in the upper half-plane; however we can not simultaneously arrange that τ is in the usual fundamental region for SL(2, ℤ), and this needs to be remembered when evaluating the Weierstrass functions below.

Of the l + 1 subgroups of ℂ/Λ of order l, the two defined over ℝ are the one generated by z = λ_{1}/l (in both cases), and in the first case, the one generated by z = λ_{2}/l, or in the second case, the one generated by z = (λ_{1} − 2λ_{2})/l. Thus z/λ_{1} is either 1/l, τ/l, or (1 − 2τ)/l. Let ℘(z; τ) denote the Weierstrass ℘-function relative to the lattice [1, τ]. Then we have

x_{k} = ℘(kzλ_{1}^{−1}; τ)λ_{1}^{−2} − \frac{1}{12}b_{4}  and  y_{k} = \frac{1}{2} (℘′(kzλ_{1}^{−1}; τ)λ_{1}^{−3} − a_{1}x_{k} − a_{3}) .

Here we have had to take account of the lattice scaling [λ_{1}, λ_{2}] = λ_{1}[1, τ], and also of the fact that (℘(z), ℘′(z)) is a point on the model of E of the form y^{2} = 4x^{3} − g_{2}x − g_{3} = 4x^{3} − (c_{4}/12)x − (c_{6}/216) rather than a standard model where the coefficient of x^{3} is 1.

We evaluate these points of order l numerically for k = 1, 2, …, (l − 1)/2, for each of the two values of z (depending on whether we are in case 1 or case 2). Substituting into Vélu’s formulae, we obtain in each case the real coefficients a′_{i} of a curve which is l-isogenous to E over ℝ. If these coefficients are close to integers we round them and check that the resulting curve over ℚ has the same conductor N as the original curve E. If not, we also test the curve with coefficients l^{i}a′_{i}.

The resulting program finds l-isogenous curves very quickly for any given prime l. We run it for all primes l in the set determined previously, applying it recursively to each new curve found until we have a set of curves closed under l-isogeny for these values of l. Since the set of primes l for which a rational l-isogeny exists is itself an isogeny invariant, once we have finished processing the first curve in the class, we will already know which primes l to use for all the remaining curves.

Some care needs to be taken with a method of computation such as this, where we use floating-point arithmetic to find integers. The series we use to compute the periods and the Weierstrass function and its derivative all converge very quickly, so that we can compute the a′_{i} to whatever precision is available, though of course in practice some rounding error is bound to arise. When we test whether a floating-point number is ‘approximately an integer’ in the program, we must make a judgement on how close is close enough. With too relaxed a test, we will find too many curves are ‘approximately integral’; usually these will fail the next hurdle, where we test the conductor, but this takes time to check (using Tate’s algorithm). On the other hand, too strict a test might mean that we miss some rational isogenies altogether, which is far more serious. In compiling the tables, there was only one case which caused trouble after the program had been finely tuned. The resulting error resulted in a curve (916B1) being erroneously listed as 3-isogenous to itself in the first (preprint) edition of the tables; this is possible only when a curve has complex multiplication, which is not the case here, though it does not often occur even in the complex multiplication case (see the remarks in the next section). Unfortunately the error was not noticed in the automatic generation of the typeset tables, and I am grateful to Elkies for spotting it.^{11} The curve E = [0, 0, 0, −1013692, 392832257] has three real points of order 2, two of which are equal to seven significant figures; the period ratio is approximately 7i. One of the curves 2-isogenous to E over ℝ has coefficients [0, 0, 0, −1013691.999999999992, 392832257.000000006], which are extremely close to those of E itself. Thus this new curve, which is not defined over ℚ, passed both our original tests (the coefficients are extremely close to integers, and the rounded coefficients are those of a curve of the right conductor, namely E itself). After this example was discovered, we inserted an extra line in the program, to print a warning whenever a supposedly isogenous curve was the original curve itself, and reran the program on all 2463 isogeny classes (which only takes a few minutes of machine time). The result was that expected, namely that 916B1 is the only curve for which this phenomenon occurs within the range of the tables^{12}. There is no example of a curve actually l-isogenous to itself with conductor less than 1000.

Our original implementation of this algorithm in Algol68 used a precision of approximately 30 significant figures for its real and complex arithmetic, which was sufficient to find all the isogenous curves up to conductor 1000. However, our implementation in C++ misses several isogenous curves when using standard double precision, with approximately 15 digits (though this runs very quickly); we need to use a multiprecision floating-point package (such as the one included in LiDIA) to obtain a satisfactory working program, though the resulting code runs very much slower. In our extended computations to conductor 5077, we have computed the isogenies independently using both a C++/LiDIA program and a PARI program, and the results agree.

When we were initially persuaded to extend the tables to include isogenous curves as well as the modular curves themselves, we were afraid that the total number of resulting curves would be rather larger than it turned out to be. On average, we found that the number of curves per isogeny class was 5113/2463, or just under 2.08. We do not know of any asymptotic analysis, or even a heuristic argument, which would predict an average number of two curves per class. However, it is dangerous to generalize from the limited amount of data which we have available. In the extended computations to conductor 5077, the ratio slowly diminishes; for all curves up to this conductor, the ratio is 31570/17583 = 1.795.

11This error also somehow survived into the first edition of this book, despite these comments in the text.

12Another example of the same type occurs for curve 1342C3, where the period ratio is approximately 9.5i.

3.9 Twists and complex multiplication

Traces of Frobenius.

If E is given by a standard minimal Weierstrass equation over ℤ, then for all primes p of good reduction the trace of Frobenius a_{p} is given by

a_{p} = 1 + p − |E(𝔽_{p})| .

If E has bad reduction at p, this same formula gives the correct value for the pth Fourier coefficient of the L-series of E.

Since in our applications we never needed to compute a_{p} for large primes p, we used a very simple method to count the number of points on E modulo p. First, for all primes p in the desired range (say 3 ≤ p ≤ 1000; p = 2 would be dealt with separately), we precompute the number n(t, p) of solutions to the congruence s^{2} ≡ t (mod p). Then we simply compute

a_{p} = p − ∑_{x=0}^{p−1} n(4x^{3} + b_{2}x^{2} + 2b_{4}x + b_{6}, p).

This was sufficient for us to compute a_{p} for all p < 1000 for all the curves in the table, which we did to compare with the corresponding Hecke eigenvalues. For large p, there are far more efficient methods, such as the baby-step giant-step method or Schoof’s algorithm (see [51]). Details of these may be found in [9]. More recently, even better algorithms have been developed, by Atkin, Elkies, Morain, Müller and others. For example, Morain and Lercier in 1995 successfully computed the number of points on the curve [0, 0, 0, 4589, 91228] over 𝔽_{p} for p = 10^{499} + 153, a prime with 500 decimal digits. This took 4200 hours of computer time.

Twists.

A twist of a curve E over ℚ is an elliptic curve defined over ℚ and isomorphic to E over \overline{ℚ} but not necessarily over ℚ itself. Thus the set of all twists of E is the set of all curves with the same j-invariant as E. These can be simply described, as follows.

First suppose that c_{4} ≠ 0 and c_{6} ≠ 0; equivalently, j ≠ 1728 and j ≠ 0 (respectively). Then the twists of E are all quadratic, in that they become isomorphic to E over a quadratic extension of ℚ. For each integer d (square-free, not 0 or 1), there is a twisted curve E ∗ d with invariants d^{2}c_{4} and d^{3}c_{6}, which is isomorphic to E over ℚ(\sqrt{d}). If E has a model of the form y^{2} = f(x) with f(x) cubic, then E ∗ d has equation dy^{2} = f(x). A minimal model for E ∗ d may be found easily by the Laska–Kraus–Connell algorithm. The conductor of E ∗ d is only divisible by primes dividing ND, where D is the discriminant of ℚ(\sqrt{d}). The simplest case is when gcd(D, N) = 1; then E ∗ d has conductor ND^{2}. More generally, if D^{2} ∤ N then E ∗ d has conductor lcm(N, D^{2}), but if D^{2} | N then the conductor may be smaller; for example, (E ∗ d) ∗ d is isomorphic to E, so has conductor N again.

Twisting commutes with isogenies, in the sense that if two curves E, F are l-isogenous then so are their twists E ∗ d, F ∗ d. If E has no complex multiplication (see below), then the structure of the isogeny class of E is a function of j(E) alone.

The trace of Frobenius of E ∗ d at a prime p not dividing N is χ(p)a_{p}, where χ is the quadratic character associated to ℚ(\sqrt{d}) and a_{p} is the trace of Frobenius of E. Thus if E is modular, attached to the newform f, then E ∗ d is also modular and attached to the twisted form f ⊗ χ, in the notation of Chapter 2.

When j = 0 (or equivalently c_{4} = 0), E has an equation of the form y^{2} = x^{3} + k with k ∈ ℤ non-zero and free of sixth powers. Such curves have complex multiplication by ℤ[(1+\sqrt{−3})/2]. Two such curves with parameters k, k′ are isomorphic over ℚ(\sqrt[6]{k/k′}).

Similarly, when j = 1728 (or equivalently c_{6} = 0), E has an equation of the form y^{2} = x^{3}+kx with k ∈ ℤ non-zero and free of fourth powers. Such curves have complex multiplication by ℤ[\sqrt{−1}]. Two such curves with parameters k, k′ are isomorphic over ℚ(\sqrt[4]{k/k′}).

Complex multiplication.

Each of the 13 imaginary quadratic orders 𝔒 of class number 1 has a rational value of j(𝔒) = j(ω_{1}/ω_{2}), where 𝔒 = ℤω_{1} + ℤω_{2}. Elliptic curves E with j(E) = j(𝔒) have complex multiplication: their ring of endomorphisms defined over ℂ is isomorphic to 𝔒. In all other cases the endomorphism ring of an elliptic curve defined over ℚ is isomorphic to ℤ, since an elliptic curve with complex multiplication by an order of class number h > 1 has a j-invariant which is not rational, but algebraic of degree h over ℚ.

We give here a table of triples (D, j, N) where j = j(𝔒) for an order of discriminant D, and N is the smallest conductor of an elliptic curve defined over ℚ with this j-invariant. All but the last three values (D = −43, −67, −163) have N < 1000 and so occur in the tables.

D −4 −16 −8 −3 −12  −27  −7 −28 −11 −19 −43 −67  −163^{j 12^{3} 66^{3} 20^{3} 0 2 · 30^{3} −3 · 160^{3} −15^{3} 255^{3} −32^{3} −96^{3} −960^{3} −5280^{3} −640320^{3}}^{N 32 32 256 27 36  27  49 49 121 361 43^{2}  67^{2}  163^{2}}

If E has complex multiplication by the order 𝔒 of discriminant D, then the twist E ∗ D is isogenous to E, though not usually isomorphic to E (over ℚ). Indeed, the only cases where E is isomorphic to E ∗ D are D = −4 and D = −16 with j(E) = 1728: the curves y^{2} = x^{3} + 16kx and y^{2} + 256kx are twists of, and isomorphic to, y^{2} = x^{3} + kx. Since curves are isogenous if and only if they have the same L-series by Faltings’s Theorem (see [22]), this implies that E has complex multiplication if and only if a_{p} = χ(p)a_{p} for all primes p, where χ is the quadratic character as above. Thus a_{p} = 0 for half the primes p, namely those for which χ(p) = −1. This gives an alternative way of recognizing a curve with complex multiplication, from its traces of Frobenius. This is particularly convenient in the case of modular curves, where we compute the a_{p} first, and will always know when a newform f, and hence the associated curve E_{f}, has complex multiplication. For, in such a case, we must have D^{2} | N and f = f ⊗ χ, which we may easily check from the tables.