Cube root computation

[johnywalkyra]

Hello,

first of all sorry for crossposting, but I could not decide which group
is more appropriate. To my question: Recently I've came across
the code in GCC standard library, which computes the cube root
of a real (floating point) number. Could anyone explain me the math
behind the computation? Here's the snippet:

---8<-----8<-----8<-----8<-----8<-----8<-----8<-----8<--

#define CBRT_2 1.2599210498948731648
#define CBRT_2_SQR 1.5874010519681994748

static const double cbrt_factors[5] =
{
1.0/CBRT_2_SQR,
1.0/CBRT_2,
1.0,
CBRT_2,
CBRT_2_SQR
};

double cbrt(double x)
{
double xm, ym, u, t2;
int xe;

xm = frexp(fabs(x), &xe);

u = (+0.354895765043919860
+ ((+1.50819193781584896
+ ((-2.11499494167371287
+ ((+2.44693122563534430
+ ((-1.83469277483613086
+ (+0.784932344976639262 -
0.145263899385486377*xm)*xm)
*xm))
*xm))
*xm))
*xm));
t2 = u*u*u;
ym = u*(t2 + 2.0*xm)/(2.0*t2 + xm)*cbrt_factors[2 + xe%3];

return ldexp(x > 0.0 ? +ym : -ym, xe/3);
}

---8<-----8<-----8<-----8<-----8<-----8<-----8<-----8<--

For those of you, who do not know frexp & ldexp functions,
here are their descriptions:

double frexp(double x, int *expptr)
The frexp function breaks down the floating-point value (x)
into a mantissa (m) and an exponent (n), such that
the absolute value of m is greater than or equal to 0.5
and less than 1.0, and x = m*2^n. The integer exponent
n is stored at the location pointed to by expptr.

double ldexp(double x, int exp)
The ldexp function returns the value of x*2^exp,
if successful.

Thanks for your concern,

John Walker jr.

[Tim Prince]

jo**********@post.cz wrote:

Recently I've came across
the code in GCC standard library, which computes the cube root
of a real (floating point) number. Could anyone explain me the math
behind the computation? Here's the snippet:

---8<-----8<-----8<-----8<-----8<-----8<-----8<-----8<--

#define CBRT_2 1.2599210498948731648
#define CBRT_2_SQR 1.5874010519681994748

static const double cbrt_factors[5] =
{
1.0/CBRT_2_SQR,
1.0/CBRT_2,
1.0,
CBRT_2,
CBRT_2_SQR
};

double cbrt(double x)
{
double xm, ym, u, t2;
int xe;

xm = frexp(fabs(x), &xe);

u = (+0.354895765043919860
+ ((+1.50819193781584896
+ ((-2.11499494167371287
+ ((+2.44693122563534430
+ ((-1.83469277483613086
+ (+0.784932344976639262 -
0.145263899385486377*xm)*xm)
*xm))
*xm))
*xm))
*xm));
t2 = u*u*u;
ym = u*(t2 + 2.0*xm)/(2.0*t2 + xm)*cbrt_factors[2 + xe%3];

return ldexp(x > 0.0 ? +ym : -ym, xe/3);
}

I don't think gcc comes with such a function; perhaps it comes with your
glibc?
This looks fairly straightforward, but unnecessarily obscure and in need
of commenting. Certainly not very original, as Googling came right up
with a 30 year old Fortran version. The argument is analyzed as xm *
pow(2,xe) * sign(x). There is a polynomial approximation, evaluated by
Horner's method, for 1/cuberoot(2) times the cube root over the range
[.5,1.), followed by a disguised Newton iteration. This result is
multiplied by cube root of pow(2,xe). It would appear that accuracy
could be improved by making the Newton iteration adjustment at the end.
For efficiency, one would want to assure that that xe/3 is not
performed twice. Note avoidance of the standard C div() function.

[Randy Poe]

jo**********@post.cz wrote:

Hello,

first of all sorry for crossposting, but I could not decide which group
is more appropriate. To my question: Recently I've came across
the code in GCC standard library, which computes the cube root
of a real (floating point) number. Could anyone explain me the math
behind the computation? Here's the snippet:

---8<-----8<-----8<-----8<-----8<-----8<-----8<-----8<--

#define CBRT_2 1.2599210498948731648
#define CBRT_2_SQR 1.5874010519681994748

static const double cbrt_factors[5] =
{
1.0/CBRT_2_SQR,
1.0/CBRT_2,
1.0,
CBRT_2,
CBRT_2_SQR
};

double cbrt(double x)
{
double xm, ym, u, t2;
int xe;

xm = frexp(fabs(x), &xe);

u = (+0.354895765043919860
+ ((+1.50819193781584896
+ ((-2.11499494167371287
+ ((+2.44693122563534430
+ ((-1.83469277483613086
+ (+0.784932344976639262 -
0.145263899385486377*xm)*xm)
*xm))
*xm))
*xm))
*xm));
As already explained, this is a polynomial approximation
for the cube root of a number in what you say is the range
[0.5,1.0]. It's a fifth-order polynomial in nested form, which
is the most efficient way to evaluate a polynomial.

When you break a number down into mantissa and exponent,
x = xm * 2^xe, the cube root is given by cbrt(xm) * 2^(xe/3.0).
However, if xe is not an integer, then 2^(xe/3.0) is going
to require some computation.

The solution is to break it down into 2^(integer part of xe/3)*
2^(fractional part of xe/3). The table handles the possible
values of 2^(1/3) and 2^(2/3).
t2 = u*u*u;
ym = u*(t2 + 2.0*xm)/(2.0*t2 + xm)*cbrt_factors[2 + xe%3];

The cbrt_factors lookup is to handle non-integer powers
of 2, as I said. I'm not sure of the purpose of the first part
of this expression. It may be a single-step Newton's
method to improve the estimate u.

- Randy