I have a small collection of utilities for working with elliptic curves. Some of these utilities (in particular, algorithms for computing rational torsion points) require solving the subproblem of finding the rational roots of a cubic polynomial with rational coefficients. The tools were written to make it easier for me to play with toy examples while studying a graduate course, so weren’t built for speed - however, I was surprised to find how difficult it was to get ahold of these roots.

I mention experiments a few times below; I have these buried in a repo somewhere, and will dig them up and tidy them one day™ when I have the time.

Here’s the specific setup. Consider \(f(X) = X^3 + aX^2 + bX + c \in \mathbb{Q}[X]\). The problem is to find all \(\alpha \in \mathbb{Q}\) such that \(f(\alpha) = 0\). We make the following, trivial, observations:

If \(\alpha\) is a rational root, then \(f(X) = (X-\alpha)g(X)\), for \(g(X) \in \mathbb{Q}[X]\) monic quadratic

\(g(X)\) has either 0 or 2 rational roots. We can efficiently derive them (and whether they exist) from the quadratic formula

So, equivalently, the problem is to find any \(\alpha \in \mathbb{Q}\) such that \(f(\alpha) = 0\).

Now, we have a simple equation for the roots of a quadratic, as mentioned. Having reasonably substantial training in field theory, I was pretty confident this would be similarly easy. We have, as a consequence of a bit of Galois theory, a pretty easy algorithm for extracting the roots of such a cubic in terms of radicals by hand, in some algebraic extension of \(\mathbb{Q}\), and so we can just use this formula to extract the roots, then just check which of these are nice rationals. This is doable by pen and paper - how hard can it be in Haskell? Absolutely awful, it turns out.

Naive attempt

In the name of sparing the reader from the required algebra, we’ll just quote the “cubic formula” (sometimes known as Cardano’s formula) here. A Galois-theoretic derivation of the general solution in an algebraic number field is quite beautiful and worth the effort to understand (in my opinion), but Wikipedia has a “bare hands” derivation of the formula here. To summarise, given \(f\) as above:

Set \(\;\Delta_0 = a^2 - 3b\), \(\;\Delta_1 = 2a^3 - 9ab + 27c\)

Set \(\;C = \sqrt[3]{\frac{\Delta_1 \pm \sqrt{\Delta_1^2 - 4\Delta_0^2}}{2}}\), where we may freely choose any cube root, and where the choice of square root is free unless one would render \(\;C=0\) - in this case, chose the root which makes \(\;C \neq 0\)

The roots of \(f\) are then given by \(\;x_k = - \frac{1}{3}\left( b + \xi^k C + \frac{\Delta_0}{\xi^k C} \right)\), for \(\;k=0, 1, 2\), where \(\;\xi\) is a 3rd root of unity

Seems easy enough - there’s this small niggle of what happens with negative square roots, but we might assume this will work itself out for cases where we have rational roots. Knocking this together, we get something like:

rationalCubicRoots :: Rational -> Rational -> Rational -> [Rational]
rationalCubicRoots a b c 
    | d0 == 0 && d1 == 0 = replicate 3 (-b/(3*a))
    | otherwise = roots where
    d0 = a^2 - 3*c :: Rational
    d1 = 2*a^3 - 9*a*b + 27*c :: Rational
    t = squareRootRational $ d1^2 - 4*d0^3 :: Maybe Rational
    c' :: Maybe Rational
    c' | d0 == 0 = cubeRootRational =<< tp
       | otherwise = cubeRootRational =<< tm where
       tp = (/2) . (d1+) <$> t
       tm = (/2) . (d1-) <$> t

    firstRoot = ((-1/3)*) . (b+) <$>
        (\x -> x + d0/x) =<< c'
    roots = maybe [] (\x -> x:rationalQuadraticRoots (a + x) (b + (a + x)*x)) rationalRoot

We quickly knock together some tests to check it does what it’s supposed to; and immediately find it fails on nearly all inputs.

The problem is, it’s not uncommon for the term \(q = d_1^2 - 4d_0^3\) to fail to be square. It’s very often negative, and this borks the whole thing. This is a well-known phenomenon, known as “casus irreducibilis” - in order to get to the rational solutions, we have to pass through intermediate terms in some complex extension of \(\mathbb{Q}\), even if the imaginary parts will eventually cancel out.

Now, this isn’t a problem in itself. It’s compounded by the fact that \(\|q\|\) may also fail to be square, so we have to work with radicals. This means symbol arithmetic. Not only that, but we then have to take a cube - while it might be plausible to keep track of the radicals involved when we’re only taking square roots and handle cancellation, once we’re taking radicals of radicals, it becomes particularly ugly.

Denesting radicals

What if we could pull out the square root? There are many well-known identities in this respect; for example, Ramanujan famously spotted:

\[\sqrt[3]{\sqrt[3]{2} - 1} = \sqrt[3]{1/3} - \sqrt[3]{2/9} + \sqrt[3]{4/9}\] \[\sqrt{\sqrt[3]{28} - \sqrt[3]{27}} = \left( \sqrt[3]{98} - \sqrt[3]{28} - 1 \right) /3\]

among others. Can we do this more generally?

The answer is yes; but not totally. A paper of Blömer gives some specific Galois-theoretic conditions for when we can denest radicals of depth two (which is our situation here), and there is some work of Landau extending this. When it is possible, it is very fiddly work on a computer; but the main takeaway is that it is not always possible.

Radicals and floating-point precision

The first thing that one might think of to circumvent this is to first compute the roots via a floating-point algorithm, and then nail down to a specific rational once we’re at high enough accuracy. Of course, \(\mathbb{Q}\) is dense in \(\mathbb{R}\), so this isn’t in itself a solution - however, we have the following:

Rational Root Theorem: Let \(\;f \in \mathbb{Z}[X]\), and let \(\;\beta \in \mathbb{Q}\) such that \(\;f(\beta) = 0\). Then, letting \(\;a, b \in \mathbb{Z}\) coprime such that \(\;\beta = \frac{a}{b}\), we have that \(\;b\) divides the leading coefficient of \(\;f\), and \(\;a\) divides the constant term of \(\;f\)

An easy corollary of this says that if we clear denominators of a rational polynomial, we can get an equivalent condition on the rational roots.

So by the rational root theorem, we know that there is only a finite set of values the denominator of a rational root of our original \(f\) can have. We can simply take an lcm of this set, giving, say, \(d\), and then we know that any rational root must be expressible as \(\frac{a}{d}\), for \(a \in \mathbb{Z}\). At this point, a reasonable line of attack is to use a floating-point root-finding algorithm with a tolerance less than \(\frac{1}{2d}\); from this we’d be able to work out \(a\) exactly.

There are broadly two ways we might try and go about this:

“Exact” methods

Use the trigonometric solutions for cubic equations, using arbitrary-precision arithmetic. In order to pull out a root, we would need to implement at least sqrt, cos, and acos. We can use MPFR, the GNU library for floating-point horrors. There are Haskell bindings in existence for this. This allows us to do floating-point computations at arbitrary precision, but we do need to establish what precision we want to work at.

So; let’s run through the analysis. The computation is described over here, and goes something like this.

Begin with the depressed cubic \(\;t^3 + pt + q = 0\). If our original polynomial is not depressed, we can obtain a depressed cubic via the substitution \(\;t \leftarrow x - \frac{a}{3}\).

Then we have cases on the number of real roots of the equation (which may be determined from the discriminant):

Case 1: 3 real roots. Then we have a root \(\;t_0 = 2 \sqrt{\frac{-p}{3}} \cos \left( \frac{1}{3} \left( \arccos{\frac{3q}{2p}\sqrt{\frac{-3}{p}}} \right)\right)\)

Let’s pause here and think about this case before we continue. Via a Taylor expansion, we can study how error in our representation of \(p\) and \(q\) will propagate. Considering just the \(\arccos\) component of this formula, and letting \(x + \delta\) being our approximation of \(x\), we have by Taylor’s theorem:

for some \(\xi\) between \(x\) and \(x + \delta\). There is a catch, though; the derivative of \(\arccos(z)\) is \(-\frac{1}{\sqrt{1- z^2}}\), which is singular within the domain we’re concerned with. This means, fixing an “acceptable error” for this step of the computation, say, \(\varepsilon\), that as the argument to \(\arccos\) gets closer and closer to \(\pm 1\), we will need to control \(\delta\) to be smaller and smaller; and that the maximal allowable value of \(\| \delta \|\) can become arbitrarily small, despite holding \(\varepsilon\) fixed.

To see this explicitly: given some \(\varepsilon\), we want \(\delta(\epsilon)>0\) such that for each \(x\), \(\| \delta(\epsilon) \cdot \frac{-1}{\sqrt{1 - (x \pm \delta(\varepsilon))^2}} \| \leqslant \varepsilon\), where we choose between \(\pm\) by taking the one closest to 1. Equivalently,

The discriminant is

and so:

Now, we have \(x \in [-1, 1]\). Fixing \(\varepsilon\) and varying \(x\), the problem becomes apparent. As \(x\) approaches \(1\), the left side of this bracket approaches \(0\) from below. As \(x\) approaches \(-1\), the right side of this bracket approaches \(0\) from above. So no matter what we set our \(\delta(\varepsilon)\) to be, there will always be a value of \(x\) in the valid range such that it isn’t small enough.

I did give this a go anyway, just to see how reasonable the bounds were. QuickCheck quickly kept finding failing counter-examples until I started running out of RAM.

Iterative methods

Search for the solution directly. In theory, this does work. In practice, all standard iteration schemes I attempted were horrifically slow to converge at this level. I did some testing with criterion to check I wasn’t imagining things, but yeah; this isn’t really an option.

Applying RRT directly

An alternative is that we just try to find the rational root in our finite set directly. The delta between this and the above method will be particularly large when the set is small, and the denominators large. This might look something like the following:

rationalCubicRoots :: Rational -> Rational -> Rational -> [Rational]
rationalCubicRoots a b 0 = 0:rationalQuadraticRoots a b
rationalCubicRoots a b c = roots where
    x3Term = foldl lcm 1 $ map denominator [a, b, c] :: Integer
    absConstTerm = abs (x3Term * numerator c) :: Integer

    generateAllFactors :: [(Prime Integer, Word)] -> [Integer]
    generateAllFactors [] = [1]
    generateAllFactors ((p,n):xs) = (take (1 + fromEnum n) . iterate (* unPrime p)) 
        =<< generateAllFactors xs

    possibleDenominators = sortOn abs $ generateAllFactors $ factorise x3Term
    factorsOfConstTerm = generateAllFactors $ factorise absConstTerm
    -- If the root is 0, it will have been picked up by the casing above
    possibleNumerators = sortOn abs $ factorsOfConstTerm ++ map negate factorsOfConstTerm

    lagrangianLimit = foldl max 1 $ map abs [a, b, c]
    possibleFractionsWithinLimit = 
        (\d -> takeWhile ((<= lagrangianLimit) . abs) $ 
            map (% d) possibleNumerators
            ) =<<
        takeWhile ((<= lagrangianLimit) . abs . (1 %)) 
        possibleDenominators

    rationalRoot = find ((== 0) . (\x -> x^3 + a*x^2 + b*x + c)) possibleFractionsWithinLimit
    roots = maybe [] (\x -> x:rationalQuadraticRoots (a + x) (b + (a + x)*x)) rationalRoot

My experiments showed this was okay; not blazing fast over test sets of thousands of cases, but at least correct. With the optimisation of computing the factors before we compute the lcm, this gives the function I’m currently using for this:

rationalCubicRoots :: Rational -> Rational -> Rational -> [Rational]
rationalCubicRoots a b 0 = 0:rationalQuadraticRoots a b
rationalCubicRoots a b c = roots where
    x3TermFactors :: Map (Prime Integer) Word
    x3TermFactors = foldl1 (Map.unionWith max) $
        map (Map.fromList . factorise . denominator) [a, b, c] 
    absConstTermFactors :: Map (Prime Integer) Word
    absConstTermFactors = Map.unionWith (+) x3TermFactors $
        Map.fromList $
        factorise $ 
        numerator c

    generateAllFactors :: [(Prime Integer, Word)] -> [Integer]
    generateAllFactors [] = [1]
    generateAllFactors ((p,n):xs) = (take (1 + fromEnum n) . iterate (* unPrime p)) 
        =<< generateAllFactors xs

    possibleDenominators = sortOn abs $ generateAllFactors $ Map.toList x3TermFactors
    factorsOfConstTerm = sortOn abs $ generateAllFactors $ Map.toList absConstTermFactors
    -- If the root is 0, it will have been picked up by the casing above
    possibleNumerators = join $ map ((^..each) . (id &&& negate)) factorsOfConstTerm

    lagrangianLimit = foldl max 1 $ map abs [a, b, c]
    possibleFractionsWithinLimit = 
        (\d -> takeWhile ((<= lagrangianLimit) . abs) $ 
            map (% d) possibleNumerators
            ) =<<
        takeWhile ((<= lagrangianLimit) . abs . (1 %)) 
        possibleDenominators

    rationalRoot = find ((== 0) . (\x -> x^3 + a*x^2 + b*x + c)) possibleFractionsWithinLimit
    roots = maybe [] (\x -> x:rationalQuadraticRoots (a + x) (b + (a + x)*x)) rationalRoot

There is potentially more performance to squeeze out by short-circuiting in cases that lend themselves nicely to Cardano — but this is more than sufficient for me, for now.