algorithm - Computing a correctly rounded / an almost correctly rounded floating-point cubic root

Question

Suppose that correctly rounded standard library functions such as found in CRlibm are available. Then how would one compute the correctly rounded cubic root of a double-precision input?

^{This question is not an “actual problem that [I] face”, to quote the FAQ. It is a little bit like homework this way. But the cubic root is a frequently found operation and one could imagine this question being an actual problem that someone faces.}

Since “best Stack Overflow questions have a bit of source code in them”, here is a bit of source code:

  y = pow(x, 1. / 3.);

The above does not compute a correctly rounded cubic root because 1/3 is not representable exactly as a double.

---

ADDITIONAL NOTES:

An article describes how to compute a floating-point cubic root, but the last iteration(s) of the recommended Newton-Raphson algorithm would have to be done to higher precision for the algorithm to compute a correctly rounded double-precision cubic root. That may be the best way to compute it, but I am still looking for a shortcut that would take advantage of existing correctly rounded standardized functions.

C99 includes a cbrt() function, but it cannot be expected to be correctly rounded or even faithful for all compilers. CRlibm's designers could have chosen to include cbrt() in the list of provided functions, but they didn't. References to implementations available in other libraries of correctly rounded math functions are welcome.

Answer 1

Given that there are lots of easily computable rational points on the curve x = y^3, I’m tempted to reduce about s^3 ~ x, with s rational and just a few bits wide. Then you have:

cbrt(x) = s * cbrt(1 + (x - s^3)/s)

Obvious thing then is to evaluate the correction term using your favorite series approximation, and compute a residual via head-tail FMA arithmetic to bump the result up or down by an ulp if needed (you won’t need the full computation most of the time, obviously).

This isn’t totally in the spirit of the question, but can definitely be made to work, and it’s very easy to prove the necessary bounds this way. Hopefully someone else can suggest something more clever (I used up my cleverness for the month already).