# Finding roots of cubic equations using iterative versus analytic techniques

The book Evaluating Derivatives states that although “the roots of cubic and even quartic equations can be expressed in terms of arithmetic operations and fractional powers, it is usually more efficient to calculate them by an iterative process”, and then goes on to describe Newton’s method. Certainly, the analytical formulae for solving cubic equations on Wikipedia look awful, but Newton’s method uses loops! Surely that’s worse.

## The iterative solution

Well, let’s see. This is implemented in my favourite programming language. We’ll start by defining Newton’s method for functions of type `f64 -> f64` in a straightforward manner, parameterised by tolerance `eps` and starting point `x0`:

``````let newton (eps: f64) (x0: f64) (f: f64 -> f64) (f': f64 -> f64) =
loop x = x0 while f64.abs(f x) > eps do x - f x / f' x``````

Then we define an evaluation function for cubic functions of the form ax³+bx²+cx+d, as well as the derivative 3ax²+2bx+c.

``````let cubic a b c d x : f64 =
a*x**3 + b*x**2 + c*x + d

let cubic' a b c _d x : f64 =
3 * a * x**2 + 2 * b * x + c``````

Now we can put these pieces together and get an iterative solver for finding a root of a cubic function:

``````let newton_cubic_root eps x0 a b c d =
newton eps x0 (cubic a b c d) (cubic' a b c d)``````

Cubic functions always have at least one real root, but may have more. This method finds one of them:

``> newton_cubic_root 1.0e-4f64 0.0f64 4.0f64 5.0f64 -6.0f64 -8.0f64``
``1.2450299590260148f64``
``> cubic 4.0f64 5.0f64 -6.0f64 -8.0f64 (newton_cubic_root 1.0e-5f64 0.0f64 4.0f64 5.0f64 -6.0f64 -8.0f64)``
``4.864979530339042e-10f64``

## The analytical solution

The following is the analytical solution. As above, it finds the first real root and ignores the others (if they exist).

``````let cubic_root a b c d : f64 =
let b = b / a
let c = c / a
let d = d / a
let q = (3*c - b*b)/9
let r = -27*d + 9*c*b - 2*(b*b*b)
let r = r / 54
let disc = q*q*q + r*r
let term1 = b/3
in if disc > 0 then
let s = r + f64.sqrt disc
let s = if s < 0 then -((-s)**(1/3)) else s**(1/3)
let t = r - f64.sqrt disc
let t = if t < 0 then -((-t)**(1/3)) else t**(1/3)
in -term1 + s + t
else if disc == 0 then
let r13 = if r < 0 then -(-r)**(1/3) else r**(1/3)
in -term1 + 2*r13
else let q = -q
let r13 = 2*f64.sqrt q
in -term1 + r13*f64.cos(q**3/3)``````

“Analytic” it may be, but given the vagaries of floating point computation, I’m not certain it’s significantly more accurate than the iterative method in practice.

``> cubic_root 4.0f64 5.0f64 -6.0f64 -8.0f64``
``1.245029959006595f64``
``> cubic 4.0f64 5.0f64 -6.0f64 -8.0f64 (cubic_root 4.0f64 5.0f64 -6.0f64 -8.0f64)``
``3.552713678800501e-15f64``

## Benchmarking yet again

Now let’s try benchmarking. Newton’s method is extremely sensitive to the starting point. I somewhat arbitrarily decide on `x0=0`. It would run substantially faster with `x0=-1` for the randomly generated coefficients we’ll use (all uniformly drawn from `[0,1]`), but that feels like cheating.

``````entry bench_newton = map4 (newton_cubic_root 0.00001 0)

entry bench_cubic_root = map4 cubic_root``````

We’ll benchmark on 100k sets of coefficients.

``````-- ==
-- entry: bench_newton bench_cubic_root
-- random input { f64 f64 f64 f64 }``````

And the results, using plain sequential execution:

``````cubic.fut:bench_newton:
f64 f64 f64 [...:       2718μs (RSD: 0.081; min:  -6%; max: +17%)

cubic.fut:bench_cubic_root:
f64 f64 f64 [...:       5574μs (RSD: 0.010; min:  -1%; max:  +2%)``````

So indeed, Newton’s method seems faster. And this is with a very naive implementation: it would be much faster (10x for `x0=-1`) if we put a little bit more effort into finding a good starting point. In contrast, while there are special cases of cubic functions that are easier to find roots for, I don’t think the general case can be solved much better than we are doing here. But who knows, maybe it can - I’m not that much of an expert here.