-- # 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](https://en.wikipedia.org/wiki/Newton%27s_method).
-- Certainly, the analytical formulae for solving cubic equations [on
-- Wikipedia](https://en.wikipedia.org/wiki/Cubic_equation#Depressed_cubic)
-- 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](https://futhark-lang.org). 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 0.0001 0.0 4.0 5.0 (-6.0) (-8.0)
-- > cubic 4.0 5.0 (-6.0) (-8.0) (newton_cubic_root 0.00001 0.0 4.0 5.0 (-6.0) (-8.0))
-- ## 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.0 5.0 (-6.0) (-8.0)
-- > cubic 4.0 5.0 (-6.0) (-8.0) (cubic_root 4.0 5.0 (-6.0) (-8.0))
-- ## 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.