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

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)

-- > 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).

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)

-- > 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.