Fifth root for doubles.

fifth_root.c

static double
fifth_root_pos_finite(double x)
{
    static const double fifth_roots_pow2[9] = { // 2^((i-4)/5)
        0.57434917749851754908974044155911542475223541259766,
        0.65975395538644709958475687017198652029037475585937,
        0.75785828325519899451023775327485054731369018554687,
        0.87055056329612412469032278750091791152954101562500,
        1.00000000000000000000000000000000000000000000000000,
        1.14869835499703509817948088311823084950447082519530,
        1.31950791077289419916951374034397304058074951171875,
        1.51571656651039798902047550654970109462738037109375,
        1.74110112659224824938064557500183582305908203125000,
    };

    double z = x;
    int e, new_e;

    if (x == 0.0)
        return x;

    // extract power of 2, reducing mantissa to [0.5,1)
    x = frexp(x, &e);

    // Approximate fifth root of number in [0.5,1), relative error below 7.78e-6
    x =   0.59493167584123773927728962007677182555198669433594 + x *
        ( 0.84132549636479991850279702703119255602359771728516 + x *
        (-0.7812833955684315156986485817469656467437744140625  + x *
        ( 0.45952167932621135193471673119347542524337768554687 + x *
        (-0.114503232341356855905623035596363479271531105041504))));

    // adjustment based on exponent%5 (covers both positive and negative exponent range)
    new_e = e/5;
    x *= fifth_roots_pow2[e - new_e*5 + 4];

    // multiply by power of 2
    x = ldexp(x, new_e);

    // refine via Newton
    x -= (x - (z / (x*x*x*x))) * 0.2;
    x -= (x - (z / (x*x*x*x))) * 0.2;

    return x;
}