#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>

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

typedef union _f32
{
    float f;
    uint32_t i;
} _f32;

float fft_fifthroot(float x)                                                    
{
    _f32 tmp;
    float t, r, y;

    tmp.f = fabsf(x);
    tmp.i   = (0x4c28cf99 - (tmp.i/5));

    printf("x0: %#.7f, S: %#.10e, ", x, 1.0f/tmp.f);
    y = tmp.f;
    t  = y * y; t *= t; t *= y;
    r = (0.2f * y * (1.0f - (x * t)));   // 1st iteration
    printf("iter0: %#.7f, r: %#.7e, x1: %#.7f, ", r, y, r+y);
    y += r;
    t  = y * y; t *= t; //t *= y;
    r  = (0.2f * (y - ((y*y) * x * t)));   // 2nd iteration
    printf("iter1: %#.7f, r: %#.7e, x2: %#.7f\n", r, y, r+y);
    y += r;
    //t  = y * y; t *= t; t *= y;
    //y += (0.2f * y * (1.0f - (x * t)));   // 3rd iteration
    return y;
}

#ifdef TEST

#include <stdio.h>

static __attribute__((always_inline)) inline void disable_denormals(void)
{
    unsigned int a = 0x9f80;
    __asm__ __volatile__ ("ldmxcsr %0" : "+m" (a));
}

int main(int argc, char **argv)
{
    float in = (float)strtod(argv[1], NULL);
    float d1, d2, d3;

    disable_denormals();
    setvbuf(stdout, NULL, _IONBF, 0);

    //d1 = fifth_rootf_pos(in);
    d2 = 1.0f/fft_fifthroot(in);
    d3 = powf(in, 0.2);

    printf("%#.7g^(1/5): %#.10f (%a) / %#.10g (%a)\n", in, d2, d2, d3, d3);
    return 0;
}

#endif