MATHEMATICS OF COMPUTATION Volume 00, Number 0, Xxxx XXXX, Pages 000–000 S 0025-5718(XX)0000-0

EFFICIENT SOLUTION OF RATIONAL CONICS

J. E. CREMONA AND D. RUSIN

Abstract. We present efficient algorithms for solving Legendre equations over ℚ (equivalently, for finding rational points on rational conics) and parametrizing all solutions. Unlike existing algorithms, no integer factorization is required, provided that the prime factors of the discriminant are known.

1. Introduction

1.1. Summary of results. In this paper we give efficient methods of finding all rational points on a rational conic C, given by non-singular homogeneous equation of degree 2:

C : f(X, Y, Z) = 0. (1)

One method for finding one rational point on C, if one exists, is the original descent method of Legendre. We show how one may easily make a significant improvement to this (reducing the number of iterations from exponential in the size of the input to linear); and also, but with more work, make an even greater improvement. This last method involves no integer factorization other than that of the discriminant of the original equation (which is in any case necessary for deciding the solubility of (1)). It is the necessity of factorizing “spurious” integers arising during the course of the computation which is the bottleneck in simpler reduction methods; our “factorization-free” method avoids this entirely.

We also describe a factorization-free method of solution based on lattice reduction; this is not original, though apparently not well known.

We present examples and timings of our implementation of both methods; these indicate that the reduction method is faster in practice than the lattice-based method. Both are linear time, given a so-called solubility certificate (defined below), and probabilistic polynomial time given only the factorization of the discriminant.

As an example of the speed which is now attainable, the solution of an equation of the form ax^{2} + by^{2} = cz^{2}, where a, b and c are 200-digit primes, takes less than 2 seconds on a modest PC. Such a problem is not feasible to solve in reasonable time with Legendre’s method (as in Maple, for example).

We also show how to parametrize all rational points on C, given one point, in the most efficient way. This is necessary for several applications, such as to 2-descent on elliptic curves, and is also used for finding a small single solution to (1).

It would be useful and interesting to extend the algorithms presented here to number fields. We say little more about this here, but refer to the paper [11] by Pohst, and Simon’s thesis [14].

The factorization-free algorithm presented here has been implemented in release 2.8 (July 2001) of the package Magma.

\begin{matrix}c \\ ◯ \end{matrix}1997 American Mathematical Society

1991 Mathematics Subject Classification. Primary 11G30, 11D41.

1.2. Background. By the Hasse or local-global principle for curves of genus 0, the curve C has rational points if and only if it has points everywhere locally. Thus, testing (1) for solubility is easy, at least in theory, and in practice no harder than factorizing the discriminant of the given equation (see Section 1 below for details).

Our first main concern will be to find efficiently one solution when solutions exist. Here and throughout we will pass freely between the geometric language of “points on curves” and the Diophantine language of “solutions to equations”. We always exclude the trivial solution (x, y, z) = (0, 0, 0), as we are really interested in projective solutions (x : y : z) ∈ ℙ^{2}(ℚ), each of which has a “primitive” representation with x, y, z ∈ ℤ and gcd(x, y, z) = 1, unique up to sign.

Secondly, we will want to find a “small” solution. Holzer’s Theorem (see below for a precise statement) asserts that a soluble equation always has solutions which are not too large in terms of the coefficients. Any given solution may be reduced, using a method of Mordell, until it satisfies Holzer’s bounds. We present an alternative reduction method, faster than Mordell’s, though the solution it gives may not be quite Holzer-reduced.

Finally, given one solution P_{0} = (x_{0}, y_{0}, z_{0}) to (1), one can write down a parametrization of all solutions of the form

X = Q_{1}(U, V), Y = Q_{2}(U, V), Z = Q_{3}(U, V) (2)

where each Q_{i}(U, V) ∈ ℤ[U, V] is a quadratic form. Geometrically, the homogeneous coordinates (U : V) parametrize the pencil of lines through P_{0}, each of which intersects the conic C in a unique second point. Our final task will be to find such a parametric solution which is as simple as possible. We will see that a parametrization exists such that the discriminants of the polynomials Q_{i}(U, V) are prescribed in terms of the coefficients of the defining polynomial f(X, Y, Z), independently of the particular basic solution P_{0} found earlier. This last point is particularly significant in certain applications.

Our approach throughout will be algorithmic, and our results will be in the form of efficient algorithms to carry out the tasks we have just outlined. We will give examples to show that our method is more efficient, and leads to better (meaning smaller) solutions than those which can be found elsewhere (for example, by using the Maple computer algebra system). The mathematics here is entirely elementary, and mostly also quite well-known, but we are not aware of a systematic treatment of such equations in the literature which is both algorithmic and concerned with the size of the parametric solutions obtained.

A slightly different problem is to parametrize all “primitive” integer solutions (x, y, z) to (1) using integer quadratics Q_{i}(U, V). Mordell showed that this is possible using a finite family of quadratic parametrizations of the form (2). Since we are interested in projective solutions we are not interested in the primitivity, and our task is therefore slightly simpler.

The application which led us to develop these methods is in higher descents on elliptic curves over ℚ, starting with a descent via 2-isogeny. See [5] for details of this. Another application of which we are aware is the determination of explicit equations for hyperelliptic curves whose Jacobians are quotients of modular Jacobians: see the theses of Wilson [19] (Oxford, 1998) and Weber [17] (Essen, 1996) for examples of these. It is remarkable that an algorithmic solution to the problem of finding all rational points on a curve of genus 0 has not yet been perfected (as remarked by Mazur in [9]), given the current activity on a wide scale concerning constructive solutions to Diophantine equations of higher genus, so it is interesting to note that efficient solutions to this simpler problem are also required for the study of curves of higher genus.

We are grateful to Denis Simon for the reference [12].

2. Single Solutions

2.1. Standard forms of equation. By elementary algebra we may complete the square in the general quadratic form f(X, Y, Z) to obtain the diagonal form

aX^{2} + bY^{2} + cZ^{2} = 0, (3)

often called Legendre’s equation. Since the equation is assumed to be nonsingular, we have abc ≠ 0. Furthermore, by simple scaling of the variables we may reduce to the case where the coefficients are integers which are (i) pairwise coprime, and (ii) square-free, so that abc is square-free. Achieving condition (i) only requires gcd computations, while (ii) requires factorization of the coefficients. We will assume throughout that this factorization is known. Such an equation (3) will be called reduced; it is unique, up to permutations of the variables and changing all the signs. Since real solubility requires that the coefficients do not all have the same sign, we also assume that a > 0, b > 0 and c < 0.

It will also sometimes be useful to put our equation into norm form

X^{2} − aZ^{2} = bY^{2}. (4)

Solving this amounts to expressing b as a norm from ℚ(\sqrt{a}), if possible. In this form we can require that a and b are both square-free integers, but not that they are coprime. Real solubility requires that a and b are not both negative. We will use this form for the first recursive solution of the equation below.

Lastly, for the applications to elliptic curves it is most convenient to use a form of the equation slightly more general than the diagonal form, which we call the semi-diagonal form:

aX^{2} + bXZ + cZ^{2} = dY^{2}. (5)

Here we will require that all the coefficients are integers with d squarefree, d(b^{2}−4ac) nonzero for nonsingularity, and gcd(a, b, c, d) = 1. In our application we also have ac ≠ 0 and so we will also assume this below.

2.2. Local Solubility Criterion and Holzer’s Theorem. The necessary and sufficient criterion for solubility of (3) is simply that it should have solutions in ℚ_{p} for all primes p and also in ℝ. This result is usually referred to as Legendre’s Theorem. For odd primes p not dividing abc there is always a local solution, so this only gives a finite number of conditions to check. Checking these conditions in practice does require us to factorize the coefficients. Suppose that (3) is reduced, so that abc is square-free. If p is odd and divides c (say), then solubility in ℚ_{p} follows from solubility modulo p (by Hensel’s Lemma), and hence from the condition that the Legendre symbol (\frac{−ab}{p}) is +1. Hence the local conditions for all odd finite primes are equivalent to the existence of solutions to the following quadratic congruences:

(6) X_{1}^{2} ≡ −bc (mod a); X_{2}^{2} ≡ −ca (mod b); X_{3}^{2} ≡ −ab (mod c).

Moreover, the number of local conditions which fails must be even (by the product formula for the Hilbert symbol), so the solubility of these congruences together with the sign condition ensuring real solubility are already sufficient to ensure global solubility, and a 2-adic condition is not needed.

Definition 2.1. A triple (k_{1}, k_{2}, k_{3}) ∈ ℤ^{3} is called a solubility certificate for (3) if it gives a solution to the congruences (6).

We summarize the local solubility criterion as follows.

Lemma 2.1. Let a, b and c be nonzero integers with abc squarefree, not all of the same sign. Then (3) has a solution if and only if a solubility certificate exists.

If a, b and c are pairwise coprime (but not necessarily square-free), then the existence of a solubility certificate is sufficient, but no longer necessary, for the existence of solutions to (3).

A proof of the last statement is implicit in the algorithms below, which guarantee to deliver a solution from a solubility certificate provided only that a, b and c are pairwise coprime and not all of the same sign. That the existence of the certificate is not necessary when the coefficients are square-free may be seen from the equation 9X^{2} − Y^{2} − Z^{2} = 0, which has the solution (1, 3, 0) but no certificate since the congruence X_{1}^{2} ≡ −1 (mod 9) has no solution.

Associated to the triple of coefficients (a, b, c) and the certificate (k_{1}, k_{2}, k_{3}) we will associate a 3-dimensional sublattice L = L(a, b, c; k_{1}, k_{2}, k_{3}) of ℤ^{3}, called the solution lattice for the certificate, as follows:

L(a, b, c; k_{1}, k_{2}, k_{3}) = {(x, y, z) ∈ ℤ^{3} |by ≡ k_{1}z (mod a),

cz ≡ k_{2}x (mod b), (7)

ax ≡ k_{3}y (mod c)}.

The index of L(a, b, c; k_{1}, k_{2}, k_{3}) in ℤ^{3} is |abc|. One easily checks that for (x, y, z) ∈ L, we have ax^{2} + by^{2} + cz^{2} ≡ 0 (mod abc). In the second and third algorithms we present below, we will construct a solution of (3) which lies in the solution lattice for any given solubility certificate.

The first algorithm we give below for solving conics itself constitutes a proof of Legendre’s theorem, since it is guaranteed to find a solution unless either a quadratic congruence fails to be soluble, or the signs of the coefficients are wrong. Indeed, Legendre’s own proof follows the same lines: see the account in Weil’s historical book [18, p. 100]. Algorithmic solutions in the literature often follow essentially the same reduction procedure as Legendre (see [8], or [15] for a recent example). As we will see, this method has two disadvantages in practice: it takes many steps, each of which involves the factorization of an integer, and the resulting solution can be very large. Our first improvement already performs better in these respects; although it does not eliminate the factorization from each step, the number of steps is reduced, the numbers to be factorized are smaller, and the resulting solution is also smaller. Then the “factorization-free” version of the reduction method eliminates the need for any factorization, given a solubility certificate, giving even greater improvement and making possible the solution of equations whose coefficients have hundreds of digits in only a few seconds.

In the famous paper [1], in which higher descents were used to study the ranks of elliptic curves of the form Y^{2} = X^{3} − DX, the authors remark [1, p. 100] that the solution of various auxiliary conics is the most time-consuming part of the descent process. We also found this to be true (despite having 30 years of factorization technology to hand) before using the new methods described here.

Now assume that (3) is soluble. Holzer’s theorem asserts that there exists an integral solution (x, y, z) with

|x| ≤ \sqrt{|bc|};  |y| ≤ \sqrt{|ac|};  |z| ≤ \sqrt{|ab|}, (8)

or equivalently,

max(|a|x^{2}, |b|y^{2}, |c|z^{2}) ≤ |abc|. (9)

Such a solution we will call “Holzer-reduced”. Holzer’s Theorem is not trivial to prove: see [3] for a recent fairly short proof, improving earlier versions by Mordell and Cassels (see section 2.4 below for more on this). In Mordell’s proof, one obtains a solution which does not necessarily satisfy Holzer’s bounds, and then reduces the solution using the following lemma from [10, Theorem 5, p.47].

Lemma 2.2. Let a, b and c be nonzero integers with abc squarefree, a > 0, b > 0 and c < 0, and let (x_{0}, y_{0}, z_{0}) be a solution of (3). If |z_{0}| > \sqrt{ab} then there exists a solution (x_{1}, y_{1}, z_{1}) with |z_{1}| < |z_{0}|.

We will give Mordell’s construction below. After a finite number of steps, we arrive at a solution (x, y, z) with |z| ≤ \sqrt{ab}, and then the inequalities on x and y follow immediately. We will also present a new method of reducing solutions which is faster than Mordell’s, but only produces a solution satisfying

max(|a|x^{2}, |b|y^{2}, |c|z^{2}) ≤ \frac{4}{3}|abc|. (10)

A similar result concerning small solutions to Legendre’s equation over totally real number fields can be found in [11].

2.3. Algorithm I: Legendre-type reduction. The first algorithm for finding one solution works with the equation in the norm form (4), where the coefficients a and b are square-free nonzero integers, not necessarily coprime. By symmetry we may assume that |a| ≤ |b|, interchanging a and b if necessary. The idea, which originates with Legendre, is to proceed by descent, reducing the problem of solving (4) to that of solving a similar equation with a smaller b coefficient. This step is repeated until |b| < |a|, after which a and b are interchanged. The base cases in which no further descent is necessary are trivially dealt with. The full procedure is as follows.

Algorithm I.

1. If |a| > |b| then swap a and b, solve the resulting equation, then swap y and z in the solution obtained.

2. If b = 1 then set (x, y, z) = (1, 1, 0) and stop.

3. If a = 1 then set (x, y, z) = (1, 0, 1) and stop.

4. If b = −1 there is no solution (since a must be −1).

5. If b = −a then set (x, y, z) = (0, 1, 1) and stop.

6. If b = a then let (x_{1}, y_{1}, z_{1}) be a solution of X_{1}^{2} + Z_{1}^{2} = aY_{1}^{2}, set (x, y, z) = (ay_{1}, x_{1}, z_{1}) and stop.

7. Let w be a solution to X^{2} ≡ a (mod b) with |w| ≤ |b|/2, and set (x_{0}, z_{0}) = (w, 1), so that x_{0}^{2} − az_{0}^{2} ≡ 0 (mod b).

8. Use lattice reduction to find a new nontrivial solution (x_{0}, z_{0}) to the congruence X^{2} − aZ^{2} ≡ 0 (mod b), with x_{0}^{2} + |a|z_{0}^{2} as small as possible.

9. Set t = (x_{0}^{2} − az_{0}^{2})/b, and write t = t_{1}t_{2}^{2} with t_{1} square-free.

10. Let (x_{1}, y_{1}, z_{1}) be a solution to X^{2} − aZ^{2} = t_{1}Y^{2}; then

(x, y, z) = (x_{0}x_{1} + az_{0}z_{1}, t_{1}t_{2}y_{1}, z_{0}x_{1} + x_{0}z_{1})

is a solution to (4): stop. By the end of step 6 we have reduced the problem to solving equations in which |b| ≥ 2, |b| > |a| and a ≠ 1 (though a = −1 is possible). The reduction step proceeds by first solving the quadratic congruence

X^{2} ≡ a (mod b)

to obtain a solution w with |w| ≤ |b|/2. The usual algorithm for this step involves factorizing b, finding a square root of a modulo each prime divisor of b, and combining them with the Chinese Remainder Theorem. All these square roots must exist if the equation passes the local solubility criterion. We then have w^{2} − a = bt, where the integer t satisfies

|t| < \frac{1}{4}|b| + 1 ≤ \frac{1}{2}|b|;

(here we use 1 ≤ |a| < |b|). The standard algorithm found in the literature (as in [15], for example) omits step 8, using the fact that this value of t is strictly less than b to obtain a descent. This procedure works perfectly well in practice provided that the initial coefficients a and b are fairly small. The size of the larger coefficient is reduced by a factor of approximately 4 at each step; the main problem with large examples is the need to factorize all the coefficients b which arise, in order to solve the associated quadratic congruences.

Our improvement consists of inserting the extra step 8 above. We have one solution (x_{0}, z_{0}) = (w, 1), to the congruence

X^{2} − aZ^{2} ≡ 0 (mod b). (11)

Using an elementary lattice reduction technique we find the solution (x_{0}, z_{0}) to this congruence which minimizes x^{2} + |a|z^{2}, and set t = (x_{0}^{2} − az_{0}^{2})/b. This will be very much smaller than the earlier value of t. Explicitly, the minimal vector has length O(b\sqrt{a}), so we see that in Step 9, t will be O(\sqrt{a}). Thus while the unimproved method only reduces the size of ab (measured in bits, say) at a rate linear in the number of steps, in the improved method the size is reduced quadratically. One expects that the number of digits in ab should be roughly halved with each iteration. We give an example in the next section.

The rest of the procedure (steps 9 and 10) is identical, with or without the lattice reduction step 8. The formula in step 10 comes from the multiplicativity of the norm from ℚ(\sqrt{a}) to ℚ:

(x_{0} + z_{0}\sqrt{a})(x_{1} + z_{1}\sqrt{a}) = (x_{0}x_{1} + az_{0}z_{1}) + (x_{0}z_{1} + z_{0}x_{1})\sqrt{a},

and hence

b(t_{1}t_{2}y_{1})^{2} = (bt_{1}t_{2}^{2})(t_{1}y_{1}^{2}) = (x_{0}^{2} − az_{0}^{2})(x_{1}^{2} − az_{1}^{2})
= (x_{0}x_{1} + az_{0}z_{1})^{2} − a(x_{0}z_{1} + z_{0}x_{1})^{2}.

Note that in Step 9, it is not really necessary to factorize t, since t_{1} need not be square-free in Step 10; but since solving the reduced equation in Step 10 will first involve factorizing t_{1} to solve the congruence X^{2} ≡ a (mod t_{1}), there is no time lost in finding this square-free decomposition immediately.

The square-free decomposition is the main time-consuming step in the algorithm, together with the solution of the subsidiary quadratic congruences in Step 7. It involves factorization of the numbers t which arise during the course of the computation, but which need not be related in any direct way to the coefficients of the original equation. We have developed a way of avoiding this factorization altogether, which will be described below in Algorithm II. Starting with a solubility certificate, either the solubility certificate at each level will determine a solubility certificate at the next level (which immediately gives the solution to the quadratic congruence we need), or alternatively a square factor of one of the coefficients will be obtained, which also leads to a reduced problem. See Section 2.5 below.

A similar idea of using 2-dimensional lattice reduction to solve a modular version of our problem was described in the paper [12].

For completeness we give the details of the lattice reduction used in Step 8. Define a positive definite quadratic form on ℤ^{2} by

(u, v) · (u′, v′) = (wu + bv)(wu′ + bv′) + |a|uu′,

so that the (square) norm of (u, v) is

∥(u, v)∥^{2} = (wu + bv)^{2} + |a|u^{2}.

Let (u_{0}, v_{0}) be the nonzero vector in ℤ^{2} which minimizes this norm. One may find (u_{0}, v_{0}) by starting with the standard basis (1, 0), (0, 1) and applying Gaussian reduction. Then set (x_{0}, z_{0}) = (u_{0}w + bv_{0}, u_{0}): we have

x_{0}^{2} − az_{0}^{2} ≡ u_{0}^{2}(w^{2} − a) ≡ 0 (mod b),

and x_{0}^{2} + |a|z_{0}^{2} = ∥(u_{0}, v_{0})∥^{2} is minimal.

2.3.1. Example. To illustrate the dramatic improvement which the lattice reduction trick (Step 8 of the algorithm) provides in a non-trivial example, we take the equation (4) with a = −113922743 and b = 310146482690273725409, which occurs in [17]. The unimproved algorithm (omitting Step 8) proceeds with 18 reduction steps and the following sets of coefficients:

(a, b) = (−113922743, 310146482690273725409)
(a, b) = (−113922743, 6322888267334211334) (a, b) = (−5941135, 690379)
(a, b) = (−113922743, 22155222796709666)  (a, b) = (690379, −5941135)
(a, b) = (−113922743, 13176519068967)  (a, b) = (690379, −436439)
(a, b) = (−113922743, 552039370818)  (a, b) = (−436439, 690379)
(a, b) = (−113922743, 10830811819)  (a, b) = (−436439, 52017)
(a, b) = (−113922743, 52527821)  (a, b) = (52017, −436439)
(a, b) = (52527821, −113922743)  (a, b) = (52017, −14)
(a, b) = (52527821, −5941135)  (a, b) = (−14, 52017)
(a, b) = (−5941135, 52527821)  (a, b) = (−14, 942)

At this stage, the congruence X^{2} + 14 ≡ 0 (mod 942) yields the solution x = 92, and luckily 92^{2} + 14 = 942t, with t = 9. As this is a square, we obtain a solution to the equation at this level. Passing back up the stack we finally obtain the following solution to the original equation:

(x : y : z) = (17096570497733995340458855914415817266660083175129

: 971656516633305795680905979479465911216

: 67668402208023840270008872724333068943397229).

By contrast, with the improved method we obtain the following much shorter sequence of coefficients:

(a, b) = (−113922743, 310146482690273725409)
(a, b) = (−113922743, 339)
(a, b) = (339, −113922743)

The last equation has solution (31006 : 1 : 1781) and two back-substitutions lead to the solution (320832774821087 : 21372 : −18438099853) of the original equation, considerably smaller than the previous solution found. Notice the dramatic reduction in the size of b at the first descent step compared with the first solution. The congruence X^{2} ≡ −113922743 (mod 310146482690273725409) has solution w = −88566846089432467791, leading to t = 25291553069336845336; but then lattice reduction finds the solution (x_{0}, z_{0}) = (824644660421, −93793135) to the congruence X^{2} ≡ −113922743Z^{2} (mod 310146482690273725409), which yields the much smaller value t = 5424 = 4^{2} · 339.

2.4. Improving the solution. The method of the preceding section will find one solution (x_{0}, y_{0}, z_{0}) to a diagonal equation (3), but this solution is not necessarily “Holzer-reduced”. It is possible to reduce the size of a solution. We present two methods for this: the first, due to Mordell in [10, Theorem 5, p.47], is guaranteed to produce a Holzer-reduced solution after a finite number of steps, but is slow in practice since the number of steps appears to be linear in the size of z_{0}. The second method is based on the quadratic parametrizations which will be introduced in Section 3. This is much faster in practice. The solution it produces is not always Holzer-reduced, as it is only guaranteed to satisfy (10) rather than (9), though in practice it usually is.

As well as applying one or other of these reduction procedures to the solution produced at the end of the recursion, it is also possible to reduce all the intermediate solutions used in the back-substitution Step 10. This can be beneficial in practice for large problems, since otherwise the exponential growth in the size of intermediate solutions can cause serious degradation of the running time owing to the need to work with very large integers.

2.4.1. Mordell’s method for reducing solutions. Mordell’s method is used in [10] to prove that Holzer-reduced solutions always exist. It is not presented there as an algorithm, but is easily turned into one. We refer to [10, Theorem 5, p.47] for the proof that this method works (as stated in Lemma 2.2 above), giving here only a sketch.

Suppose that we have a primitive solution (x_{0}, y_{0}, z_{0}) to the equation (3), where abc is squarefree, a and b are positive and c is negative, with |z_{0}| > \sqrt{ab}. Since c is squarefree, we have gcd(x_{0}, y_{0}) = 1.

If c is even, set k = c/2, solve k = uy_{0} −vx_{0} for u and v, and let w be the nearest integer to −(aux_{0} + bvy_{0})/(cz_{0}). Then the equations

x = \frac{1}{k} (x_{0}(au^{2} + bv^{2} + cw^{2}) − 2u(aux_{0} + bvy_{0} + cwz_{0}))
(12)  y = \frac{1}{k} (y_{0}(au^{2} + bv^{2} + cw^{2}) − 2v(aux_{0} + bvy_{0} + cwz_{0}))
z = \frac{1}{k} (z_{0}(au^{2} + bv^{2} + cw^{2}) − 2w(aux_{0} + bvy_{0} + cwz_{0}))

define integers x, y, z which also satisfy (3), and which satisfy 0 < |z| < |z_{0}|. (These follow easily from the identity (aux_{0} + bvy_{0} + cwz_{0})^{2} + ab(uy_{0} − vx_{0})^{2} = −kcz_{0}z, together with the inequalities |aux_{0} + bvy_{0} + cwz_{0}| ≤ \frac{1}{2}|cz_{0}| and ab < z_{0}^{2}.)

If c is odd, solve c = uy_{0} − vx_{0} for u and v. Now let w be the nearest integer to −(aux_{0} + bvy_{0})/(cz_{0}) which has the same parity as au + bv. Then (12) with k = 2c again defines an integral solution to (3), with 0 < |z| < |z_{0}|.

If the new z is still too big, we apply this again; after a finite number of steps the Holzer bounds (8) will be satisfied.

When we apply this in practice, we may either just apply it once, at the end, or alternatively we may apply it to each solution in the recursive stack before backsubstituting at Step 10.

2.4.2. Reducing solutions via quadratic parametrization. Starting from a primitive solution (x_{0}, y_{0}, z_{0}) to a diagonal equation (3) we apply the method of Section 3 below to obtain three parametrizing quadratic polynomials Q_{i}(U, V) with respective discriminants −4bc, −4ac and −4ab. As shown in the proof of Corollary 3.2 below, after applying Gaussian reduction to whichever one of the ±Q_{i} is a positive definite quadratic form, we obtain a parametrization whose leading coefficients give a solution (x_{1}, y_{1}, z_{1}) to (3) satisfying the “almost-Holzer” bounds (21). See Section 3 below for details.

When we apply this in practice, we have found that in most cases, the new solution obtained does in fact satisfy the Holzer bounds. Exceptional cases arise when the root of the reduced positive definite quadratic, which lies in the usual fundamental domain for SL(2, ℤ) in the upper half-plane, has imaginary part less than one.

For example, the equation X^{2} + 3Y^{2} = 91Z^{2} has the solution (19, 1, 2) which is not Holzer-reduced since x^{2} = 361 > 273 = 3 · 91 and z^{2} = 4 > 3 = 1 · 3. The parametrizing quadratics are X = 19U^{2} − 16UV − 11V^{2}, Y = U^{2} − 20UV + 9V^{2} and Z = 2(U^{2} − UV + V^{2}), with minimal discriminants 1092 = −4 · 3 · (−91), 364 = −4 · 1 · (−91) and −12 = −4 · 1 · 3 respectively. The latter is positivedefinite and reduced, with root (−1 + i\sqrt{3})/2, and this method cannot reduce the solution further. However, applying Holzer’s method once to (19, 1, 2) gives the Holzer-reduced solution (4, 5, 1). The corresponding parametrizing quadratics are X = 4U^{2} + 30UV − 12V^{2}, Y = 5U^{2} − 8UV − 15V^{2} and Z = U^{2} + 3V^{2}.

Note also that the leading coefficients of the reduced parametrizing quadratics need not necessarily be coprime, so the solution (x_{1}, y_{1}, z_{1}) may not be primitive; if not, we obviously obtain a further reduction by cancelling the common factor.

2.4.3. Example. Continuing our earlier example, the solution (320832774821087 : 21372 : −18438099853) is not Holzer-reduced, but applying Mordell reduction once yields the Holzer-reduced solution (30106379962113 : 7913 : 12747947692).

The much larger solution produced by the unimproved algorithm requires 27 steps of Mordell reduction to obtain the Holzer-reduced solution (47464775475069 : 3131 : 2629196804).

Using the quadratic parametrization method to reduce the solutions, we obtain the new solutions (7523107023591 : 7244 : 11931641701) (starting from the smaller original solution) and (70647575606369 : 5679 : 6632499416) (starting from the larger). These solutions are both Holzer-reduced.

Note that, as illustrated by these examples, there is nothing at all canonical in the solutions obtained, even amongst those which satisfy the Holzer bounds. The solution obtained will depend on all the choices of modular square roots made along the way, each such choice being equally valid, and leading to a distinct solution. In fact, one remarkable feature of Holzer’s theorem (apparent from its proof) is that it not only guarantees one reduced solution, but one in each class of modular solutions modulo abc, the number of which is around 2^{k} where k is the number of odd prime factors of abc.

2.5. Algorithm II: factorization-free reduction method. The preceding algorithm is adequate for solving equations where the coefficients are of “reasonable” size: reasonable in the sense that numbers of this size may be factorized quickly. But for larger problems, the time taken for intermediate factorizations make it impractical. For example, if we take the coefficients in (3) to be primes of around 100 digits (chosen so that (3) is soluble) then in the second step of the recursion one is likely to have to factorize a random integer with between 90 and 100 digits.

To avoid this, we have developed an alternative method which is quite similar in theory but avoids all factorization. The idea is that given a solubility certificate (k_{1}, k_{2}, k_{3}) for the diagonal equation (3) with coefficients (a, b, c), we can use it to construct a new solubility certificate for a reduced problem with smaller coefficients, without any (further) factorization, together with a linear transformation mapping solutions of the reduced problem to solutions of the original. While this idea is simple in principle, complications arise in practice, since at the general stage we cannot assume that the various triples of integers which arise as coefficients are square-free.

One starts with the data (a, b, c; k_{1}, k_{2}, k_{3}) and recursively reduces this to a similar set of data which is smaller, in a suitably defined sense, until one reaches an easy case where the solution can be written down immediately or found directly. The coefficients a, b and c are assumed to be pairwise coprime, but not necessarily square-free. The essential idea is that whenever certain numbers which would be coprime in the square-free case are encountered and found not to be coprime, the common factor can be used to reduce the problem by giving a non-trivial square factor of one of a, b or c, which may then be divided out.

At each stage we take care to provide a solution which lies in the lattice defined by the current solubility certificate, so that at the end the solution obtained lies in the lattice defined by the original certificate.

The next lemma will deal with the base cases under this scheme.

Lemma 2.3. If two of the coefficients of equation (3) are ±1, then a solution in the solution lattice may be found from a solubility certificate in time O(log |abc|).

Proof. By symmetry we may assume that ab = ±1. By a changing the sign of all three coefficients, and also of the solubility certificate (in order to keep the solution lattice unchanged) we may also assume that a = 1. The certificate consists of an integer k = k_{3} satisfying k^{2} ≡ −b = ±1 (mod c).

If a = 1 and b = −1 then we have the trivial solution (1, 1, 0), but we must find a solution satisfying x ≡ ky (mod c) where k is a fixed square root of +1 modulo c. If k ≡ ±1 (mod c) we simply use (1, ±1, 0). Otherwise, let c^{+} = gcd(k − 1, c) and c^{−} = gcd(k + 1, c) and set z = c^{+}c^{−}/c. One may check that in all cases z = ±1 or z = ±2; this is straightforward when c is square-free but needs a little care when 4 | c. Now the required solution is (x, y, z) = (\frac{1}{2}(c^{−} − zc^{+}), \frac{1}{2}(c^{−} + zc^{+}), z) which duly satisfies x^{2} − y^{2} + cz^{2} = 0 and x ≡ ky (mod c).

Now assume that a = b = 1, and so c < 0. Let x + yi = gcd(k + i, c) in the Euclidean ring ℤ[i] of Gaussian integers. Then it is easy to see (by considering the prime factorization of c) that x^{2} + y^{2} = |c|, so that (x, y, 1) is a solution to (3), and satisfies x ≡ ky (mod c) as required. The gcd may be computed in O(log |c|) steps by [13].

The general reduction step will start with a triple of coefficients (a, b, c), pairwise coprime but not necessarily square-free, defining an equation (3), together with a solubility certificate (k_{1}, k_{2}, k_{3}). We then construct a new equation

a′(X′)^{2} + b′(Y′)^{2} + c′(Z′)^{2} = 0, (3)′

with smaller coefficients (a′, b′, c′), a new solubility certificate (k′_{1}, k′_{2}, k′_{3}), and a linear map T from the new solution lattice L′ to L, mapping solutions to (3)′ to solutions to (3).

During the main reduction step, it can happen that we find a non-trivial square factor u^{2} of one of the coefficients. The following trivial lemma may then be used to reduce the problem; however, it is not possible to do so in such a way as to preserve the solution lattice. For this reason we do not in fact use this lemma in our implementation.

Lemma 2.4. Let (k_{1}, k_{2}, k_{3}) be a solubility certificate for (3), where the coefficients (a, b, c) are pairwise coprime. Suppose that u^{2} | a for some integer u > 1. Let u′ be an inverse of u modulo bc. Then (k_{1}, u′k_{2}, u′k_{3}) is a solubility certificate for the equation with coefficients (a/u^{2}, b, c). Also, if (x′, y′, z′) is a solution to the latter equation, then (x′, uy′, uz′) is a solution to the original equation.

Proof. Trivial.

As an example to show that this lemma cannot be strengthened to respect the solution lattices in all cases, take the equation p^{2}X^{2} + Y^{2} = Z^{2} with certificate (k_{1}, 0, 0) satisfying k_{1}^{2} ≡ 1 (mod p^{2}), and take u = p. Then u′ = p also, and the reduced equation is (X′)^{2} + (Y′)^{2} = (Z′)^{2} with the same certificate (k_{1}, 0, 0). The new solution lattice is the whole of ℤ^{3}; given a solution (x′, y′, z′) satisfying (x′)^{2} + (y′)^{2} = (z′)^{2} and the vacuous condition y′ ≡ k_{1}z′ (mod 1), the solution to the original equation is (x, y, z) = (x′, py′, pz′). This satisfies y − k_{1}z ≡ 0 (mod p), but not 0 (mod p^{2}).

The following lemma is crucial: it shows that we may find partial factorizations of any finite set of nonzero integers which approximate a full square-free decomposition, using only the operations of gcd and exact integer division.

Lemma 2.5. Let a_{i} for 1 ≤ i ≤ n be nonzero integers. There exist integers b_{i} for 1 ≤ i ≤ n, and pairwise coprime integers c_{I} indexed by the nonempty subsets I ⊆ {1, 2, …, n}, such that for 1 ≤ i ≤ n we have

a_{i} = b_{i}^{2} ∏_{I, i∈I} c_{I}. (13)

Moreover, these integers may be computed from the a_{i} using only the operations of gcd and exact integer division, in O(∑ log(a_{i})) steps.

Proof. We initialize by setting each b_{i} = 1, c_{{i}} = a_{i}, and the other c_{I} = 1. Then (13) is satisfied, but the coprimality conditions may not hold.

If gcd(c_{I}, c_{J}) = d > 1 for two subsets I and J, we divide c_{I} and c_{J} by d, multiply c_{I+J} by d (where I +J is the symmetric difference (I ∪J)−(I ∩J)), and multiply b_{i} by d for all i ∈ I ∩J. This preserves the relations (13) while decreasing the product of all the c_{I} by a factor of d. Hence, after a finite number of steps we achieve the conditions stated.

Remark. Of course, without the last sentence of its statement, the lemma would be trivial, using the prime factorizations of the a_{i}, and we could even require that the c_{I} should be square-free. A useful trick in practice is to use a small amount of trial division at the start: to ensure that the c_{I} are not divisible by the square of any prime p ≤ p_{0}, say, we may (for each such p in turn) divide out the largest even power of p from c_{{i}} and adjust b_{i} accordingly.

We will apply this lemma with n = 3 below.

Now we come to the main reduction step, which constructs a new reduced equation together with a solubility certificate, and an appropriate linear transformation.

Proposition 2.6. Given data (a, b, c; k_{1}, k_{2}, k_{3}) with (k_{1}, k_{2}, k_{3}) a solubility certificate for (3) and |bc| > 1, with (a, b, c) pairwise coprime and not all of the same sign. There is an algorithm, requiring no factorization, which either finds a non-trivial square factor of one of a, b, or c, or constructs a smaller set of data (a′, b′, c′; k′_{1}, k′_{2}, k′_{3}) such that (k′_{1}, k′_{2}, k′_{3}) is a solubility certificate for the equation (3)′, together with a linear transformation from solutions of this equation to solutions of (3).

Remark. When we use the algorithm described below, the situation where we fail to construct a reduced equation only arises when one of the coefficients is not squarefree. At the top level we will insist that the coefficients are square-free, so this cannot happen there; this is reasonable, since the criterion for solubility given in Lemma 2.1 requires square-free coefficients. At lower levels, if we identify a nontrivial square factor of one of the coefficients, then rather than deal with this using Lemma 2.4, we instead pass the factor we have found back to the level above.

Proof. We will subdivide the proof into a number of steps.

Step 1: Preliminaries. Let w ≡ c^{−1}k_{1} ≡ −bk_{1}^{−1} (mod a). Consider the sublattice of ℤ^{2} defined by the congruence y ≡ wx (mod a), with ℤ-basis (1, w), (0, a), together with the weighted Euclidean norm ||(x, y)|| = |b|x^{2} + |c|y^{2}. Let (w_{1}, w_{2}) be a minimal nonzero vector in this lattice (obtained by Gaussian reduction). Then

k_{1}w_{1} ≡ cw_{2}, k_{1}w_{2} ≡ −bw_{1} (mod a),

so that

bw_{1}^{2} + cw_{2}^{2} = at (14)

with t a small integer. Explicitly, 0 < ||(w_{1}, w_{2})|| ≤ \frac{4}{π} |a|\sqrt{|bc|}, so we have |t| ≤ \frac{4}{π} \sqrt{|bc|}. Note that by minimality of the vector (w_{1}, w_{2}), we know that gcd(w_{1}, w_{2}) has no prime factors p which do not divide a, for then p^{2} | t and we could divide out a factor p^{2} from (14). In particular, gcd(w_{1}, w_{2}) is coprime to bc.

To ease notation we use the abbreviations (u, v) = gcd(u, v), and write u ⊥ v to mean gcd(u, v) = 1.

Step 2: The reduced coefficients. Using Lemma 2.5, applied to the three integers bc, a, t, and using the fact that a ⊥ bc, we may write

bc = α^{2}b′c′,  a = β^{2}n_{1}n_{3},  t = γ^{2}n_{2}n_{3}c′

where the integers n_{1}, n_{2}, n_{3}, b′ and c′ are pairwise coprime. If |α| > 1 then either u = (α, b) or u = (α, c) gives a non-trivial square divisor u^{2} of b or c respectively, and we may stop. So we may assume α = 1. Similarly, we may assume β = 1, since otherwise we have a non-trivial square factor of a.

Hence the above equations simplify to

bc = b′c′,  a = n_{1}n_{3},  t = γ^{2}a′c′

where a′ = n_{2}n_{3}. Both triples (a′, b′, c′) and (n_{1}, n_{2}, n_{3}) are pairwise coprime. Moreover, a′, b′ and c′ cannot all have the same sign, since then t > 0 and bc > 0; but then (14) would imply that a, b, c all had the same sign.

Define d_{1} = (c, c′) and d_{2} = (b, c′). Then c′ = ±d_{1}d_{2} with d_{1} | c, d_{2} | b, and (d_{1}, d_{2}) = 1. Adjust the sign of d_{1} or d_{2} if necessary so that c′ = d_{1}d_{2}.

Remark. When we call this Proposition recursively with the reduced coefficients, it will return either a solution to the reduced equation or a non-trivial square factor of one of a′, b′ or c′. In the former case we transform the solution using Step 5 below, returning the result to the level above (or stopping if we are already at the top level). If we obtain a factor f^{2} of a′, then we divide a′ by f^{2} and multiply γ by f and repeat the process at this level. Similarly if we obtain a square factor of c′. Finally, if we obtain a non-trivial square factor f^{2} of b′, then at least one of gcd(f, b) and gcd(f, c) will be greater than 1, and can be passed back as a non-trivial square factor of b or c respectively. This final possibility cannot happen at the top level, where the coefficients are square-free. Clearly the number of these “back-tracking” steps will be finite.

Step 3: Refinements. In this step we show that various coprimality and divisibility conditions can be assumed between these variables, since otherwise non-trivial square factors of a, b or c are found. This will enable us to construct the new solubility certificate in Step 4.

• d_{i} | w_{i} for i = 1, 2:

For d_{i} | c′ | t, hence (14) implies that d_{i} | w_{i}^{2}. Let e = (d_{i}, w_{i}), and write d_{i} = eu and w_{i} = ev with u ⊥ v. Then d_{i} | w_{i}^{2} ⟹ u | ev^{2} ⟹ u | e ⟹ u^{2} | d_{i}. This gives a non-trivial square factor of either b or c, unless u = 1, so we may assume that u = 1 and deduce that d_{i} | w_{i} for i = 1, 2.

Now we may divide by c′ = d_{1}d_{2} in (14) to obtain

d_{1} \frac{b}{d_{2}} (\frac{w_{1}}{d_{1}} )^{2} + d_{2} \frac{c}{d_{1}} (\frac{w_{2}}{d_{2}} )^{2} = aa′γ^{2}. (15)

• (b, γ) = (c, γ) = 1:

For let d = (b, γ). Then d ⊥ c, so (14) implies that d | w_{2}^{2}. As before, this implies d | w_{2} (else we obtain a non-trivial square factor of b). Then (14) implies d^{2} | bw_{1}^{2}. But (b, w_{1}, w_{2}) = 1 by the remarks made in Step 1, so d^{2} | b and we have a square factor of b unless d = 1. Similarly, (c, γ) = 1, else we have a square factor of c.

• (d_{2}, w_{1}) = (d_{1}, w_{2}) = 1:

For let d = (d_{2}, w_{1}). Then d ⊥ d_{1}, so d | (w_{1}/d_{1}). Now (15) implies d | aa′γ^{2}. But d is coprime to each of a, γ and a′, since d | d_{2} | b | b′c′, so d = 1. Similarly, (d_{1}, w_{2}) = 1.

• (w_{1}, n_{2}) = (w_{2}, n_{2}) = 1:

For let d = (w_{1}, n_{2}). Then d | cw_{2}^{2} from (14), but d ⊥ c since n_{2} ⊥ bc, so d | w_{2}^{2}. Now a prime divisor p of d would divide both w_{1} and w_{2}, but not divide a = n_{1}n_{3}, contradicting the observation made below (14). Hence d = 1. The proof that (w_{2}, n_{2}) = 1 is similar.

Step 4: The new certificate. Next we define the new solubility certificate (k′_{1}, k′_{2}, k′_{3}) for the equation with coefficients a′ = n_{2}n_{3}, b′ = (c/d_{1})(b/d_{2}), c′ = d_{2}d_{1} as follows.

(16)  k′_{1} = \begin{cases}−bw_{1}w_{2}^{−1} & (mod n_{2}) \\ −k_{1} & (mod n_{3}) \end{cases}

(17)  k′_{2} = \begin{cases}(aγ)^{−1}k_{3}w_{1} & (mod c/d_{1}) \\ (aγ)^{−1}k_{2}w_{2} & (mod b/d_{2}) \end{cases}

(18)  k′_{3} = \begin{cases}k_{2}a′γw_{1}^{−1} & (mod d_{2}) \\ −k_{3}a′γw_{2}^{−1} & (mod d_{1}) \end{cases}

This uniquely defines the k′_{i} modulo a′, b′, c′ respectively by the Chinese Remainder Theorem, since (n_{2}, n_{3}) = (c/d_{1}, b/d_{2}) = (d_{1}, d_{2}) = 1. By Step 3, all the modular inverses in these formulae exist.

To check that we do indeed have a solubility certificate is now straightforward; each of the required quadratic congruences is proved in two steps using the factorization of the relevant modulus. Note that at present the signs of the k′_{i} are immaterial; but they will be important in Step 5.

• (k′_{1})^{2} + b′c′ ≡ 0 (mod a′): First, modulo n_{2} we have

(k′_{1})^{2} + b′c′ ≡ (bw_{1}w_{2}^{−1})^{2} + bc ≡ bw_{2}^{−2}(bw_{1}^{2} + cw_{2}^{2}) ≡ 0,

since n_{2}|(bw_{1}^{2} + cw_{2}^{2}); note that it also follows that k′_{1} ≡ cw_{2}w_{1}^{−1} (mod n_{2}). Modulo n_{3} we have

(k′_{1})^{2} + b′c′ ≡ k_{1}^{2} + bc ≡ 0,

since n_{3}|a.

• (k′_{2})^{2} + a′c′ ≡ 0 (mod b′): First, modulo c/d_{1} we have

(k′_{2})^{2} + a′c′ ≡ (aγ)^{−2}k_{3}^{2}w_{1}^{2} + a′c′
≡ (aγ)^{−2}(k_{3}^{2}w_{1}^{2} + a^{2}t)
≡ a^{−1}γ^{−2}(−bw_{1}^{2} + at) ≡ 0,

while modulo b/d_{2} we have

(k′_{2})^{2} + a′c′ ≡ (aγ)^{−2}k_{2}^{2}w_{2}^{2} + a′c′
≡ a^{−1}γ^{−2}(−cw_{2}^{2} + at) ≡ 0.

• (k′_{3})^{2} + a′b′ ≡ 0 (mod c′): Observe that the equation (14) implies that

(c/d_{1})w_{2}^{2} ≡ aa′d_{2}γ^{2} (mod d_{1})

on dividing by d_{1} and then reducing modulo d_{1}. Now, modulo d_{1} we have

d_{2}w_{2}^{2}((k′_{3})^{2} + a′b′) ≡ d_{2}(k_{3}a′γ)^{2} + a′b′d_{2}w_{2}^{2} ≡ −d_{2}ab(a′γ)^{2} + a′b′d_{2}w_{2}^{2}
≡ −a′b(c/d_{1})w_{2}^{2} + a′b′d_{2}w_{2}^{2} ≡ a′w_{2}^{2}(b′d_{2} − bc/d_{1})
≡ 0,

since bc/d_{1} = b′c′/d_{1} = b′d_{2}. This implies that (k′_{3})^{2}+a′b′ ≡ 0 (mod d_{1}) since (d_{2}w_{2}^{2}, d_{1}) = 1. A similar calculation shows that (k′_{3})^{2} + a′b′ ≡ 0 (mod d_{2}).

Step 5: The linear transformation. Recall that we defined in (7) a lattice L = L(a, b, c; k_{1}, k_{2}, k_{3}) associated to equation (3) and its solubility certificate. Let L′ = L(a′, b′, c′; k′_{1}, k′_{2}, k′_{3}) be the similarly defined lattice for the reduced data. We now define a linear transformation T : L′ → L, which maps solutions to the new equation into solutions to the original.

Given (x′, y′, z′) ∈ L′, set T (x′, y′, z′) = (x, y, z) where

x = −γn_{3}x′;
(19)  y = \frac{1}{n_{2}} (\frac{c}{d_{1}} \frac{w_{2}}{d_{2}} y′ + w_{1}z′);
z = \frac{1}{n_{2}} (\frac{b}{d_{2}} \frac{w_{1}}{d_{1}} y′ − w_{2}z′) .

Assuming for the moment that y, z ∈ ℤ and that T (L′) ⊆ L, direct calculation shows that

n_{2}^{2} (ax^{2} + by^{2} + cz^{2}) = aa′γ^{2} (a′(x′)^{2} + b′(y′)^{2} + c′(z′)^{2}) .

Hence T maps solutions of the equation a′(x′)^{2} + b′(y′)^{2} + c′(z′)^{2} = 0 to solutions of the equation ax^{2} + by^{2} + cz^{2} = 0. Nontrivial solutions are mapped to nontrivial solutions, since T has nonzero determinant; specifically, another direct calculation shows that |det T | = (γn_{3})^{3}n_{1}/n_{2} ≠ 0.

Note that in general T (L′) ≠ L, since

[L : T (L′)] = [ℤ^{3} : L′][L′ : T (L′)]/[ℤ^{3} : L]
= a′b′c′(γn_{3})^{3}n_{1}/(n_{2}abc)
= γn_{3}^{3}.

If t were square-free and coprime to a, we would have γ = n_{3} = 1. In this situation, which is usually the case in practice, we do have T (L′) = L.

It remains to show that (19) does define a well-defined map from L′ to L.

• y, z ∈ ℤ: Since n_{2}y ∈ ℤ and n_{2} ⊥ c′, it suffices to show that c′(n_{2}y) ≡ 0 (mod n_{2}). But modulo n_{2} we have

c′(n_{2}y) ≡ cw_{2}y′ + c′w_{1}z′ ≡ w_{1}(k′_{1}y′ + c′z′) ≡ 0,

since k′_{1}y′ + c′z′ ≡ 0 (mod a′) (from the definition of L′) and n_{2}|a′, and we have used cw_{2} ≡ k′_{1}w_{1} (mod n_{2}). The verification that z ∈ ℤ is similar.

• by ≡ k_{1}z (mod a): using n_{2} ⊥ a and c′ ⊥ a, we compute modulo a:

n_{2}c′(k_{1}z − by) ≡ k_{1}(bw_{1}y′ − c′w_{2}z′) − b(cw_{2}y′ + c′w_{1}z′)
≡ by′(k_{1}w_{1} − cw_{2}) − c′z′(k_{1}w_{2} + bw_{1}) ≡ 0.

• cz ≡ k_{2}x (mod b): Here it suffices to work modulo d_{2} and b/d_{2} separately,

since they are coprime, and we may multiply by n_{2} since n_{2} ⊥ b. Modulo d_{2}, we have

n_{2}(k_{2}x − cz) ≡ −a′k_{2}γx′ − c (\frac{b}{d_{2}} \frac{w_{1}}{d_{1}} y′ − w_{2}z′)
≡ −a′k_{2}γx′ − b′w_{1}y′ (since d_{2}|w_{2})

≡ −a′k_{2}γx′ + w_{1}k′_{3}x′ (since k′_{3}x′ ≡ −b′y′ (mod c′) and d_{2}|c′)

≡ 0 (since k′_{3}w_{1} ≡ k_{2}a′γ).

Modulo b/d_{2}:

n_{2}(k_{2}x − cz) ≡ −a′k_{2}γx′ + cw_{2}z′

≡ k′_{2}z′k_{2}γ + cw_{2}z′ (since −a′x′ ≡ k′_{2}z′ (mod b′) and (b/d_{2})|b′)

≡ a^{−1}k_{2}^{2}w_{2}z′ + cw_{2}z′
≡ a^{−1}w_{2}z′(k_{2}^{2} + ac) ≡ 0.

• ax ≡ k_{3}y (mod c): We work modulo d_{1} and c/d_{1} separately.

Modulo d_{1}, we have ax ≡ k_{3}y ⟺ k_{3}x ≡ −by, and

n_{2}(k_{3}x + by) ≡ −k_{3}a′γx′ + b (\frac{c}{d_{1}} \frac{w_{2}}{d_{2}} y′)
≡ w_{2}(k′_{3}x′ + b′y′) ≡ 0.

Modulo c/d_{1}:

n_{2}(k_{3}y − ax) ≡ k_{3}(w_{1}z′) + aγa′x′
≡ aγ(k′_{2}z′ + a′x′) ≡ 0.

When implementing this method, we first factorize the coefficients of the given diagonal equation, removing square factors and common factors of the coefficients. Then we compute a solubility certificate, returning failure if none exists. A recursive procedure based on Lemma 2.3 and Proposition 2.6 is used to find a solution in the solution lattice defined by the certificate. Finally, the solution is reduced in size using the algorithms of Section 2.4.2 and (if necessary) Section 2.4.1.

This completes the description of the factorization-free reduction method.

2.6. Algorithm III: lattice methods. Both the preceding algorithms use 2dimensional lattice reduction. One can also use the 3-dimensional lattice L = L(a, b, c; k_{1}, k_{2}, k_{3}) (defined in (7)) directly as follows. As already observed, for (x, y, z) ∈ L we have f(x, y, z) = ax^{2} + by^{2} + cz^{2} ≡ 0 (mod abc). Moreover, Minkowski’s Theorem implies that L contains a nonzero vector (x, y, z) satisfying Holzer’s bounds (8). This implies that

−|abc| < f(x, y, z) < 2|abc|,

so that either f(x, y, z) = 0 or f(x, y, z) = |abc|. In the former case, we have a Holzer-reduced solution to (3), but in the latter case we do not have a solution. Various ways of fixing this problem have been proposed, by Mordell, Cassels and more recently by Cochrane and Mitchell in [3]. They impose extra 2-adic conditions to define a sublattice L′ of index 2 in L, such that points (x, y, z) ∈ L′ satisfy f(x, y, z) ≡ 0 (mod 2abc), and apply a theorem of Gauss to assert the existence of a point (x, y, z) ∈ L′ with |f(x, y, z)| < 2|abc|, giving a solution. The case a = b = 1 requires special treatment in the proof, as in Lemma 2.3.

In order to turn this into an algorithm for solving equations in practice, one needs methods of finding short vectors in 3-dimensional lattices, since the shortest vector in L′ certainly gives a solution. In most cases, the first vector in an LLLreduced basis of L gives a solution, and the following lemma^{1} says that one does not have to look much further.

Lemma 2.7. Let b_{1}, b_{2}, b_{3} be an LLL-reduced basis of a 3-dimensional lattice L. Then the shortest vector of L has the form n_{1}b_{1} + n_{2}b_{2} + n_{3}b_{3} where each n_{i} ∈ {−1, 0, 1}.

Remark. Since ±v have the same length, this leaves us with 13 nonzero vectors to check to find the shortest vector, given an LLL-reduced basis.

Proof. Let b_{i}^{∗} for i = 1, 2, 3 be the orthogonalized basis vectors in ℝ^{3}, so that

b_{1} = b_{1}^{∗},
b_{2} = b_{2}^{∗} + μ_{21}b_{1}^{∗},
b_{3} = b_{3}^{∗} + μ_{31}b_{1}^{∗} + μ_{32}b_{2}^{∗},

Since the b_{i} are LLL-reduced, we have |μ_{ij}| ≤ 1/2 for 1 ≤ j < i ≤ 3, and

|b_{i}^{∗}|^{2} ≥ (\frac{3}{4} − μ_{i,i−1}^{2}) |b_{i−1}^{∗}|^{2} ≥ \frac{1}{2}|b_{i−1}^{∗}|^{2}

for i = 2, 3. Hence

|b_{1}^{∗}|^{2} ≤ 2|b_{2}^{∗}|^{2} ≤ 4|b_{3}^{∗}|^{2}.

Let x = α_{1}b_{1} + α_{2}b_{2} + α_{3}b_{3} ∈ L with α_{i} ∈ ℤ. Then

|x|^{2} = (α_{1} + μ_{21}α_{2} + μ_{31}α_{3})^{2}|b_{1}^{∗}|^{2} + (α_{2} + μ_{32}α_{3})^{2}|b_{2}^{∗}|^{2} + α_{3}^{2}|b_{3}^{∗}|^{2}.

Suppose |x| < |b_{1}|. We will show that this forces |α_{i}| ≤ 1 for i = 1, 2, 3. First of all, from α_{3}^{2}|b_{3}^{∗}|^{2} ≤ |x|^{2} < |b_{1}^{∗}|^{2} ≤ 4|b_{3}^{∗}|^{2} we have α_{3}^{2} < 4, so |α_{3}| ≤ 1.

Case 1: α_{3} = 0. Then α_{2}^{2}|b_{2}^{∗}|^{2} ≤ |x|^{2} < |b_{1}^{∗}|^{2} ≤ 2|b_{2}^{∗}|^{2} implies α_{2}^{2} < 2, so |α_{2}| ≤ 1. If α_{2} = 0, then x = α_{1}b_{1}, so α_{1} = 0. Otherwise, α_{2} = ±1, giving |b_{1}^{∗}|^{2} > |x|^{2} ≥ (α_{1} + μ_{21}α_{2})^{2}|b_{1}^{∗}|^{2}, which implies (α_{1} ± μ_{21})^{2} < 1, so |α_{1}| ≤ 1 since |μ_{21}| < \frac{1}{2} and α_{1} ∈ ℤ.

Case 2: α_{3} = ±1. Now

|b_{1}^{∗}|^{2} > |x|^{2} ≥ (α_{2} + α_{3}μ_{32})^{2}|b_{2}^{∗}|^{2} + |b_{3}^{∗}|^{2}

≥ \frac{1}{2}(α_{2} ± μ_{32})^{2}|b_{1}^{∗}|^{2} + \frac{1}{4}|b_{1}^{∗}|^{2},

so (α_{2} ± μ_{32})^{2} ≤ \frac{3}{2}, giving |α_{2}| ≤ 1.

We now give a recipe for bases of the lattices L and L′, in terms of the solubility certificate (k_{1}, k_{2}, k_{3}).

1shown to us by Robin Chapman

Basis for L: Using the fact that a, b, c are pairwise coprime, solve the following for u, v, a′ and b′:

ub + vc = 1; aa′ + bcb′ = 1.

Now set

α ≡ b′ck_{1} (mod a);  β ≡ ua′bk_{3} (mod bc);  γ ≡ va′ck_{2} (mod bc).

The following vectors give a basis for L:

v_{1} = (bc, 0, 0);  v_{2} = (aβ, a, 0);  v_{3} = (αβ + γ, α, 1);

for one easily checks that these vectors all satisfy the defining congruences for L, and they evidently generate a lattice of index |abc|.

Basis for L′: One easily checks that the map ϵ : L → ℤ/2ℤ given by (x, y, z) ↦ f(x, y, z)/(abc) (mod 2) is an additive homomorphism. It is surjective, since the images of (bc, 0, 0), (0, ac, 0) and (0, 0, ab) (which all lie in L) are bc, ac and ab (mod 2) respectively, and at least one of these is odd. Hence L′ = {v ∈ L | f(v) ≡ 0 (mod 2abc)} is a sublattice of L of index 2. [Again, we are grateful to R.J.Chapman for this observation.]

Let v_{i} be a basis vector of L (from the above list) such that ϵ(v_{i}) = 1 (mod 2)). Define

w_{j} = \begin{cases}2v_{i} & if j = i; \\ v_{j} − v_{i} & if j ≠ i and ϵ(v_{j}) = 1; \\ v_{j} & if j ≠ i and ϵ(v_{j}) = 0. \end{cases}

Then w_{1}, w_{2}, w_{3} is a basis for L′.

Now use a standard integer LLL-algorithm, such as in [4], to find an LLL-reduced basis b_{1}, b_{2}, b_{3} of L′ with respect to the norm ||(x, y, z)||^{2} = |a|x^{2} + |b|y^{2} + |c|z^{2}. Then for at least one of the 13 nonzero vectors v = n_{1}b_{1} + n_{2}b_{2} + n_{3}b_{3} (up to sign) we have f(v) = 0 by Lemma 2.7.

It would also be possible to use the algorithm of Vallée (see [16]) for finding the shortest vector in the 3-dimensional lattice L′. We have not implemented this.

2.7. Other methods. Finally we mention that there are two methods for solving Legendre’s equation due to Gauss: see [7, Arts. 294, 295]. These both involve the theory of reduction of ternary quadratic forms: specifically, in both solutions one constructs an indefinite ternary form of determinant −1 and reduces it to the form x^{2} + 2yz using a suitable unimodular substitution. While Gauss does give an algorithm for this reduction in [7, Arts. 272, 274], it does not seem to be very efficient in practice. Without a fast method of carrying out such a reduction, Gauss’s methods of solving Legendre’s equation are much slower than the method we presented above.

3. Parametric Solutions

Now we have one solution (x_{0}, y_{0}, z_{0}) to our equation (1), and we wish to parametrize all solutions. Our starting point is a classical method (see [10]), which was also used in [6] and may also be found in the book [15].

3.1. The Diagonal Equation. First assume that our equation is in diagonal form (3) with abc square-free. Assuming that z_{0} ≠ 0 by symmetry, one sets X = x_{0}W + U, Y = y_{0}W + V, z = z_{0}W and eliminates W to obtain the following parametric solution:

x = Q_{1}(U, V) = ax_{0}U^{2} + 2by_{0}UV − bx_{0}V^{2},
(20)  y = Q_{2}(U, V) = −ay_{0}U^{2} + 2ax_{0}UV + by_{0}V^{2},
z = Q_{3}(U, V) = az_{0}U^{2} + bz_{0}V^{2}.

These quadratics have the following discriminants:

disc(Q_{1}) = −4bcz_{0}^{2};  disc(Q_{2}) = −4acz_{0}^{2};  disc(Q_{3}) = −4abz_{0}^{2}.

Also, the 3 × 3 matrix of coefficients of the Q_{i} (which is used in the application in [6]) has determinant −4abcz_{0}^{3}. We claim that the powers of z_{0} which appear here are entirely superfluous and may be removed. This is hardly surprising, since we made an arbitrary choice of the variable Z at the start. But it is significant, since in many of the applications, such as the one in [6] and our own in 2-descent on elliptic curves, it is crucial to keep these quantities as small as possible, and to avoid introducing spurious prime factors. Our result is as follows.

Proposition 3.1. Let a, b, c be nonzero integers with abc squarefree, such that the equation aX^{2}+bY^{2}+cZ^{2} = 0 has a (nontrivial) solution. Then the set of all rational solutions may be parametrized in the form (2) where each Q_{i}(U, V) ∈ ℤ[U, V] is quadratic, with discriminants

disc(Q_{1}) = −4bc,  disc(Q_{2}) = −4ac,  disc(Q_{3}) = −4ab,

and the determinant of the coefficient matrix of the Q_{i} is 4abc. Moreover, these discriminants cannot be further reduced.

Proof. We start with the parametrization given by (20) in terms of a primitive solution (x_{0}, y_{0}, z_{0}) with z_{0} ≠ 0. It is sufficient to find an integer e such that Q_{i}(U + eV/z_{0}, V/z_{0}) has integer coefficients for i = 1, 2, 3, since this change of variables clearly reduces the discriminants of each Q_{i} by a factor of z_{0}^{2} as required, and the coefficient determinant by a factor z_{0}^{3}.

Since a is squarefree, gcd(y_{0}, z_{0}) = 1, so we can find an integer e satisfying

ey_{0} ≡ x_{0} (mod z_{0}^{2}).

From ax_{0}^{2} + by_{0}^{2} + cz_{0}^{2} = 0 it easily follows that ae^{2} ≡ −b (mod z_{0}^{2}), and then

eax_{0} + by_{0} ≡ 0 (mod z_{0}^{2}),
e^{2}ax_{0} + 2eby_{0} − bx_{0} ≡ 0 (mod z_{0}^{2}).

Now

Q_{1}(U + eV/z_{0}, V/z_{0}) = ax_{0}U^{2} + 2\frac{eax_{0} + by_{0}}{z_{0}} UV + \frac{e^{2}ax_{0} + 2eby_{0} − bx_{0}}{z_{0}^{2}} V^{2},

which has integral coefficients^{2}. A similar check shows that the coefficients of Q_{i}(U + eV/z_{0}, V/z_{0}) are also integral for i = 2 and i = 3.

For the last statement, observe that since abc is squarefree, the only square dividing all the discriminants −4ab, −4ac, −4bc is 4. Now −ab, −ac, and −bc cannot all be discriminants: none is a multiple of 4, and they cannot all be congruent to 1 (mod 4) since their product is −(abc)^{2}.

2In fact, R. Buchholz has observed that it is sufficient for e to satisfy ey_{0} ≡ x_{0} (mod z_{0}); this may produce smaller coefficients at this stage, but the reduction given in Corollary 3.2 below makes this redundant.

Corollary 3.2. With the notation as in Proposition 3.1, there exist values (u_{0}, v_{0}) of the parameters (U, V) such that gcd(u_{0}, v_{0}) = 1 and if we set x_{1} = Q_{1}(u_{0}, v_{0}), y_{1} = Q_{2}(u_{0}, v_{0}), z_{1} = Q_{2}(u_{0}, v_{0}) then (x_{1} : y_{1} : z_{1}) is a solution of (3) satisfying the “almost-Holzer” bounds

|x_{1}| ≤ \sqrt{4|bc|/3};  |y_{1}| ≤ \sqrt{4|ac|/3};  |z_{1}| ≤ \sqrt{4|ab|/3}. (21)

Proof. Assume, without loss of generality, that a > 0, b > 0 and c < 0. Then Q_{3}(U, V) is definite (and even positive definite if we take z_{0} > 0, as we may). Applying standard Gaussian reduction to Q_{3}, we find a unimodular substitution of the parameters (U, V), say U = αU′ + βV′, V = γU′ + δV′ with αδ − βγ = ±1, so that the transformed quadratic Q_{3}^{∗}(U′, V′) = Q_{3}(U, V) has leading coefficient z_{1} = Q_{3}^{∗}(1, 0) = Q_{3}(α, γ) satisfying z_{1}^{2} ≤ | disc(Q_{3})/3| = 4ab/3. Applying the same transformation to Q_{1}(U, V) and Q_{2}(U, V) we obtain new parametrizing quadratics Q_{i}^{∗} satisfying aQ_{1}^{∗}(U, V)^{2}+bQ_{2}^{∗}(U, V)^{2}+cQ_{3}^{∗}(U, V)^{2} = 0 and the same discriminants as the Q_{i}. Substituting (U, V) = (1, 0) we obtain a new solution x_{1} = Q_{1}^{∗}(1, 0), y_{1} = Q_{2}^{∗}(1, 0), z_{1} = Q_{3}^{∗}(1, 0), with z_{1}^{2} ≤ 4ab/3. Finally, ax_{1}^{2} ≤ ax_{1}^{2} + by_{1}^{2} = cz_{1}^{2} ≤ 4|abc|/3 so that x_{1}^{2} ≤ 4|bc|/3, and similarly y_{1}^{2} ≤ 4|ac|/3. This proves the result, with (u_{0}, v_{0}) = (α, γ).

3.2. Example. We apply the method of the previous section to the equation

X^{2} + 113922743 Z^{2} = 310146482690273725409 Y^{2}

treated earlier, starting with the primitive and Holzer-reduced solution (x, y, z) = (70647575606369, 5679, 6632499416). We obtain the parametrization

X = 70647575606369 U^{2} − 272768472153240 UV − 236838674874023V^{2},
Y = 5679 U^{2} − 536 UV + 20073V^{2},
Z = 6632499416 U^{2} + 24254293278UV − 24587834368V^{2}.

These have discriminants 4 · 113922743 · 310146482690273725409, −4 · 113922743 and 4 · 310146482690273725409, as expected. While the size of the coefficients in this parametrization may seem large (up to 15 digits), recall that the coefficients of the original equation have 9 and 21 digits. By comparison, the parametrization given in [17], obtained using Maple, involves coefficients all of which have between 25 and 35 digits; and more seriously, the discriminants of the quadratics given there are k^{2} times the ones given above, where k = 2^{5} · 3 · 59 · 67 · 79 · 149 · 1993 · 7187 · 45757 · 16215770450329.

3.3. The Semi-Diagonal Equation. For convenience for our elliptic curve applications, we give an alternative form of the Proposition 3.1 suited to the semidiagonal form.

Proposition 3.3. Let a, b, c, d be integers with acd(b^{2}−4ac) ≠ 0 and d squarefree, such that the equation aX^{2} + bXZ + cZ^{2} = dY^{2} has a (nontrivial) solution. Then the set of all rational solutions may be parametrized in the form (2) where each Q_{i}(U, V) ∈ ℤ[U, V] is quadratic, with discriminants

disc(Q_{1}) = 4cd,  disc(Q_{2}) = b^{2} − 4ac,  disc(Q_{3}) = 4ad.

Proof. Rather than change variables and apply Proposition 3.1 it is simpler to start from scratch with a primitive solution (x_{0}, y_{0}, z_{0}) to (5). Set Δ = b^{2} − 4ac.

We first suppose that y_{0} ≠ 0, which will certainly be the case unless Δ is a square. One parametrization is given by

X = Q_{1}(U, V) = x_{0}U^{2} + 2(bx_{0} + 2cz_{0})UV + x_{0}ΔV^{2},
(22)  Y = Q_{2}(U, V) = y_{0}U^{2}  − y_{0}ΔV^{2},
Z = Q_{3}(U, V) = z_{0}U^{2} − 2(bz_{0} + 2ax_{0})UV + z_{0}ΔV^{2},

with disc(Q_{1}) = 16cdy_{0}^{2}, disc(Q_{2}) = 4y_{0}^{2}Δ, and disc(Q_{3}) = 16ady_{0}^{2}. These must now be divided by (2y_{0})^{2}.

The argument is slightly complicated by the fact that we cannot assume that either gcd(x_{0}, y_{0}) = 1 or gcd(z_{0}, y_{0}) = 1. But gcd(x_{0}, z_{0}) = 1 since d is squarefree, so without loss of generality we may assume that x_{0} is odd, and we may factorize y_{0} = y_{1}y_{2} with gcd(2y_{1}, y_{2}) = gcd(2y_{1}, x_{0}) = gcd(y_{2}, z_{0}) = 1. Hence by the Chinese Remainder Theorem we may find e satisfying

(23)  ex_{0} ≡ −(2cz_{0} + bx_{0}) (mod 4y_{1}^{2}),
(24)  ez_{0} ≡ (2ax_{0} + bz_{0}) (mod y_{2}^{2}).

Explicitly, if sx_{0} + tz_{0} = 1 then we may set e = t(2ax_{0} + bz_{0}) − s(2cz_{0} + bx_{0}) (mod 4y_{0}^{2}). In particular, we may compute e in practice without having to determine the factorization y_{0} = y_{1}y_{2}.

Using the fact that ax_{0}^{2}+bx_{0}z_{0}+cz_{0}^{2} ≡ 0 (mod y_{0}^{2}), simple calculations show that (23) also holds modulo y_{2}^{2} and that (24) also holds modulo 4y_{1}^{2}; hence both hold modulo 4y_{0}^{2}. Also, squaring (23) and using ax_{0}^{2} + bx_{0}z_{0} + cz_{0}^{2} ≡ 0 (mod y_{0}^{2}) we find that e^{2} ≡ Δ (mod 4y_{1}^{2}), and (24) similarly implies that the same congruence holds modulo y_{2}^{2}, so we have e^{2} ≡ Δ (mod 4y_{0}^{2}). Now a trivial calculation shows that the quadratics Q_{i}(U + eV/(2y_{0}), V/(2y_{0})) have integer coefficients and the properties stated.

The case where y_{0} = 0 may easily be handled: this can only happen when Δ is a square, say Δ = δ^{2}. Note that δ ≡ b (mod 2). We start with the parametrization

x = Q_{1}(U, V) = \frac{1}{2}(ad(δ − b)U^{2} + (δ + b)V^{2}),
y = Q_{2}(U, V) = aδUV,
z = Q_{3}(U, V) = a^{2}dU^{2} − aV^{2},

with discriminants 4a^{2}cd, a^{2}Δ and 4a^{3}d respectively. Write a = a_{1}a_{2} where a_{1} = gcd(a, (δ + b)/2). Then a_{2} divides (δ − b)/2, and a simple calculation shows that the quadratics (1/a_{1})Q_{i}(U/a_{2}, V) have the desired properties.

4. Timings

We have implemented the algorithms described in Section 2 using C++ together with the LiDIA library (version 1.4) for multiprecision integer arithmetic and factorization routines. Modular square roots were computed using the implementation in LiDIA of Shanks’s RESSOL algorithm in order to find certificates. The integer LLL algorithm was implemented by us following [4].

As sample problems we considered only diagonal equations aX^{2} +bY^{2} +cZ^{2} = 0 where a, b and −c are primes chosen so that the equation is soluble. The advantage of using prime coefficients is that the factorization-free solution methods then have to do no factorization at all (other than verifying that the coefficients are prime, which was done using a standard pseudo-primality test).

We precomputed sets of test data as follows, for

k ∈ {5, 10, 15, 20, 25, 50, 75, 100, 125, 150, 175, 200, 500, 1000}.

For each k, let a be the smallest prime above 10^{k}, and b the next prime after a; then for each of the next primes p after b such that the equation aX^{2} + bY^{2} = pZ^{2} is soluble, we store the triple (a, b, c) with c = −p. The corresponding data set for each k will be denoted S_{k}. We precomputed datasets containing 100 triples for k ≤ 200, five triples for k = 500 and just one for k = 1000. This last one uses coefficients a = 10^{1000} + 453, b = 10^{1000} + 1357, and c = −(10^{1000} + 2713). We remark that computing these data sets took longer than solving the corresponding equations using our algorithms.

Each of the algorithms was then used to compute (reduced) solutions to each of the equations in each data set. For k > 20 it was not practical to use Lagrange reduction (Algorithm I), either with (LAG+R) or without (LAG) the lattice reduction improvement described above. This was partly because of the excessively long time this would have taken, but also because a bug in LiDIA’s MPQS factorization routine meant that integers of this size could often not be factored reliably, so that the timings obtained on repeated runs were very inconsistent. Since the coefficients used are prime, no factorization at all was needed for either the factorization-free reduction method (FFR) or the LLL-based methods.

The results are as follows, given in seconds, based on a DEC alpha EV6. Recall that each entry for k ≤ 200 gives the time taken to solve 100 different problems of size around 10^{k} for data set S_{k}, the datasets for k = 500 and k = 1000 having size 5 and 1 respectively.

k  LAG LAG+R  FFR  LLL
5  35.243s 4.612s 0.407s  0.362s
10 169.764s 11.869s 0.765s  0.737s
15  18.554s 1.139s  1.185s
20  31.978s 2.537s  2.629s
25  2.763s  2.982s
50  7.168s  8.839s
75  13.073s 17.819s
100  21.920s 34.871s
125  30.611s 52.856s
150  40.603s 74.219s
175  57.991s 109.221s
200  73.597s 147.364s
500  32.372s 69.576s
1000  34.031s 79.320s

We may draw the conclusion that methods which do not require factorization at intermediate stages are much faster than those which do. Of the factorization-free methods, LLL and the factorization-free reduction methods are of comparable speed for small and medium-sized problems, but for larger problems the reduction method starts to gain, being twice as fast for the problems with 200 digit coefficients.

The above comparative timings between the FFR and LLL methods are somewhat misleading, however. By using equations with prime coefficients we have avoided all factorization in computing the solutions, but both the FFR and LLL methods start by computing the solubility certificate (k_{1}, k_{2}, k_{3}) for each equation, and this computation takes a substantial proportion of the total time. To investigate further, we isolated the time for this step, which involves the computation of three square roots modulo primes of size 10^{k} for each equation, and found that it takes almost all the time of the FFR method, and about half the time of the LLL method. In the following table, the first column of timings gives the times for just computing the certificates; the next two columns give the time to find the solutions, given the precomputed certificates using both FFR and LLL methods.

k Certificate FFR LLL
5  0.237 0.181 0.175
10  0.470 0.257 0.329
15  0.735 0.327 0.522
20  1.981 0.404 0.766
25  2.093 0.484 1.025
50  5.801 0.914 3.283
75  10.809 1.411 7.290
100  18.343 2.055 16.372
125  26.485 2.746 25.966
150  35.092 3.475 38.445
175  50.607 4.633 58.593
200  63.534 5.883 83.435
500  30.003 1.888 39.575
1000  30.248 1.983 48.766

It is now apparent that our FFR method is very much faster than LLL in finding the solution from the certificate, by a factor of about 15 in the largest example. We give some more details of this last computation: 10 levels of recursion were needed; at each depth except one, the value of the variable γ is 1, the exceptional value being 5. This agrees with our expectation that γ = 1 in most cases. With the FFR method, the solution (x, y, z) produced initially had content gcd(x, y, z) = 6, with no cancelling of common factors during the recursion (in order to stay on the appropriate lattice). This small content shows the efficiency of the formulae used to map the solutions back from lower levels. After cancelling this common factor, the non-reduced solution has integers x, y, z each of 1004 digits, with “Holzer measure” max{x^{2}/|bc|, y^{2}/|ac|, z^{2}/|ab|} = 5.7 · 10^{7}. After reduction (using the quadratic parametrization) we obtain the Holzer-reduced solution with integers of 1000 digits each and Holzer measure 0.54. The LLL method produces a solution which has Holzer measure 0.47, and again 1000 digits for each of x, y and z.

In our implementation of Lemma 2.5 we use the technique mentioned in the remark after that Lemma above, to ensure that the factors c_{I} are not divisible by p^{2} for primes p < 20. For all the examples, this was sufficient to avoid ever having to backtrack, since (without this adjustment) the only square factors of the coefficients which were discovered at lower levels were products of the primes ≤ 11. A small time saving was achieved in this way.

Analyzing these two algorithms further may throw some light on this marked difference in their running times for large problems. In both cases we start by constructing a 3-dimensional lattice L in which the solution will lie. With the LLL method, we repeatedly find new bases for this same lattice, while with the FFR method we construct a new lattice at each step, and the only lattice reduction we do is on 2-dimensional projections of these. These successive lattices have decreasing index in ℤ^{3}, and their bases have smaller and smaller integer coordinates, so we would expect the computations to be faster as we go deeper into the recursion.

References

1. B. J. Birch and H. P. F. Swinnerton-Dyer, Notes on elliptic curves II, J. Reine Angew. Math. 218 (1965), 79–108.

2. J. W. S. Cassels, Lectures on elliptic curves, London Mathematical Society Student Texts, no. 24, Cambridge University Press, 1991.

3. T. Cochrane and P. Mitchell, Small solutions of the Legendre equation, Journal of Number Theory 70 (1998), 62–66.

4. H. Cohen, A course in computational algebraic number theory (third corrected printing), Graduate Texts in Mathematics, no. 138, Springer-Verlag, 1996.

5. J. E. Cremona,  Higher descents on elliptic curves,  preprint:  see

http://www.maths.nott.ac.uk/personal/jec/papers/d2.ps.

6. I. Gaál, A. Pethő, and M. Pohst, Simultaneous representation of integers by a pair of ternary quadratic forms–with an application to index form equations in quartic number fields, Journal of Number Theory 57 (1996), 90–104.

7. C. F. Gauss, Disquisitiones arithmeticae, Springer-Verlag, 1986.

8. K. Ireland and M. Rosen, A classical introduction to modern number theory, Graduate Texts in Mathematics, no. 84, Springer-Verlag, 1982.

9. B. Mazur, On the passage from local to global in number theory, Bull. Amer. Math. Soc. (N.S.) 29 (1993), no. 1, 14–50.

10. L. J. Mordell, Diophantine equations, Pure and Applied Mathematics, no. 30, Academic Press, 1969.

11. M. Pohst, On Legendre’s equation over number fields, Publ. Math. (Debrecen) 56 (2000), 535–546.

12. J. M. Pollard and C. P. Schnorr, An efficient solution of the congruence X^{2} + kY^{2} ≡ m (mod n), IEEE Transactions on Information Theory 33 (1987), no. 5, 702–709.

13. H. Rolletschek, On the number of divisions of the euclidean algorithm applied to gaussian integers, J. Symb, Comput. 2 (1986), 261–291.

14. D. Simon, Équations dans les corps de nombres et discriminants minimaux, Ph.D. thesis, Université Bordeaux I, 1998.

15. N. P. Smart, The algorithmic resolution of Diophantine equations: a computational cookbook, London Mathematical Society Lecture Notes Series, no. 117, Cambridge University Press, 1998.

16. B. Vallée, Algorithmique dans les réseaux de petite dimension: un point de vue affine sur la récherche des minima, Séminaire de Théorie des nombres de Bordeaux (1985–1986), no. 13.

17. H.-J. Weber, Algorithmische Konstruktion hyperelliptischer Kurven mit kryptographischer Relevanz und einem Endomorphismenring echt größer als Z, Ph.D. thesis, Institut für Experimentelle Mathematik, University of Essen, 1997.

18. A. Weil, Number theory: an approach through history from Hammurapi to Legendre, Birkhäuser, 1984.

19. J. Wilson, Curves of genus 2 with real multiplication by a square root of 5, Ph.D. thesis, University of Oxford, 1998.

School of Mathematical Sciences, University of Nottingham, University Park, Nottingham NG7 2RD, UK.

E-mail address: John.Cremona@nottingham.ac.uk

Department of Mathematical Sciences, Northern Illinois University, DeKalb, IL 60115, USA.

E-mail address: rusin@math.niu.edu