#include <stdio.h>
#include <math.h>
double          cbrt(double c)
{
    int             n;
    int             exponent;
    int             sign;
    double          x;
    double          x2;
    double          x6;
    double          c2;
    double          x3c;
    double          x3c16;
    double          tenx6;
    double          temp;
    double          ans0;
    double          ans1;

    sign = (c < 0.0 ? 1 : 0);
    if (sign)
        c = -c;
    c = frexp(c, &exponent);
    x = ((((1.3584464340920900529734e-1 * c
            - 6.3986917220457538402318e-1) * c
           + 1.2875551670318751538055e0) * c
          - 1.4897083391357284957891e0) * c
         + 1.3304961236013647092521e0) * c
        + 3.7568280825958912391243e-1;

    x2 = x * x;
    x6 = x2 * x2 * x2;
    c2 = c * c;
    x3c = x2 * x * c;
    x3c16 = x3c * 16.;
    tenx6 = 10. * x6;
    temp = x3c16 + tenx6 + c2;
    ans0 = x * (4.0 * (x6 + c2) + 19.0 * x3c) / (temp);
    ans1 = 1. / x2 * c * (temp) / (x3c16 + 10. * c2 + x6);
    x = (ans0 + ans1) * 0.5;
    printf("primitive guess for x is %f\n", x);
    n = (exponent % 3);
    exponent -= n;
    if (n < 0) {
        n += 3;
        exponent -= 3;
    }
    if (n == 1)
        x *= 1.25992104989487316476721; /* cube root of 2 */
    else if (n == 2)
        x *= 1.5874010519681994747517;  /* cube root of 4 */
    x = ldexp(x, exponent / 3);
    return (sign ? -x : x);
}
int             main()
{
    double          c;
    double          ans;
    c = 27;
    ans = pow(c, 1. / 3.0);
    printf("Answer using pow  = %.20f\n", ans);
    ans = cbrt(c);
    printf("Answer using cbrt = %.20f\n", ans);
    return 0;
}