# 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`

:

```
def 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*.

```
def cubic a b c d x : f64 =
a*x**3 + b*x**2 + c*x + d
def 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:

```
def 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 0.0001 0.0 4.0 5.0 (-6.0) (-8.0)`

`1.2450299590260145f64`

`> cubic 4.0 5.0 (-6.0) (-8.0) (newton_cubic_root 0.00001 0.0 4.0 5.0 (-6.0) (-8.0))`

`4.86492623963386e-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).

```
def 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.0 5.0 (-6.0) (-8.0)`

`1.245029959006595f64`

`> cubic 4.0 5.0 (-6.0) (-8.0) (cubic_root 4.0 5.0 (-6.0) (-8.0))`

`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 { [100000]f64 [100000]f64 [100000]f64 [100000]f64 }
```

And the results, using plain sequential execution:

```
cubic.fut:bench_newton:
[100000]f64 [100000]f64 [100000]f64 [...: 2718μs (RSD: 0.081; min: -6%; max: +17%)
cubic.fut:bench_cubic_root:
[100000]f64 [100000]f64 [100000]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.