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