Cube root computation

Discussion in 'C Programming' started by johnywalkyra@post.cz, Jan 16, 2006.

  1. Guest

    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.
     
    , Jan 16, 2006
    #1
    1. Advertising

  2. Tim Prince Guest

    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.
     
    Tim Prince, Jan 16, 2006
    #2
    1. Advertising

  3. Randy Poe Guest

    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
     
    Randy Poe, Jan 16, 2006
    #3
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. mrwoopey
    Replies:
    3
    Views:
    2,320
    pipodyer
    May 9, 2006
  2. Joel Leong
    Replies:
    3
    Views:
    3,202
    Scott Allen
    Apr 19, 2005
  3. sarah28
    Replies:
    3
    Views:
    1,035
    sarah28
    Jun 26, 2003
  4. Replies:
    4
    Views:
    1,042
  5. Zangief Ief

    Math cube root

    Zangief Ief, Jul 11, 2009, in forum: Ruby
    Replies:
    12
    Views:
    210
    Sergey Avseyev
    May 14, 2011
Loading...

Share This Page