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

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

/* Constants for cbrtf(x) -- single-precision cube-root of x. */
static const unsigned
B1 = 1419915941; /* B1 = (169+1/3-0.0661629)*2**23 */

static inline float rcpss(float in)
{
    float out;
    __asm__ __volatile__("rcpss %1, %0\n" :"=x"(out) :"x"(in));
    return out;
}

static inline float rcp(float in)
{
    float out;
    __asm__ __volatile__("rcpss %1, %0\n" :"=x"(out) :"x"(in));
    return ((out+out) - out*out*in);
}

float cbrtf_estimate(float x)
{
    _f32 tmp;
    tmp.f = fabsf(x); /* x <- |x| */
    /* rough cbrt to 5 bits */
    tmp.i = (B1-tmp.i/3);
    //tmp.i &= 0x7fffff00;
    return tmp.f;
}

float fft_cbrtf(float x)
{
    _f32 tmp;
    float r, t;
    uint32_t sign = (x < 0.0f);

    tmp.f = fabsf(x); /* x <- |x| */

    if(!tmp.i || tmp.i>=0x7f800000) return(x); /* cbrt(0,NaN,INF) is itself */

    /* rough cbrt to 5 bits */
    tmp.i = (B1-tmp.i/3);
    //tmp.i &= 0x7fffff00;
    t = tmp.f;

    /* new cbrt to 23 bits */
    //r = t*t;
    //t = 0.333333333f * (2.0f*(t+t) - r*x*r);
    //tmp.f = t; tmp.i &= 0x7fffff00; t = tmp.f; 
    //r = t*t;
    //t = 0.333333333f * (2.0f*(t+t) - r*x*r);
    r = t*t;
    t = 0x1.555555p-02f * (2.0f*(t+t) - r*x*r);
    //r = t*t*t;
    //t = t + (0x1.555556p-02f * (t - t*x*r));
    r = t*t*t;
    t = t + (0x1.555555p-02f * (t - t*x*r));
    //t = 0x1.555555p-02f * (4.0f*t - t*x*r);

    /* convert 1/cbrtf(x) to cbrtf(x) */
    //t *= t;
    //t *= fabsf(x);

    /* retore the sign bit */
    if (sign) t = -t;
    return(t);
}

#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <sys/time.h>

static float rel_err_diff(float x1, float x2)
{
  if(x1 == 0.0f && x2 == 0.0f) {
    return 0.0f;
  } else if (x1 == 0.0f) {
    return fabsf(x2);
  } else if(x1 <= 1e38f && x2 <= 1e38f && (x1 + x2 != 0.0f)) {
    //return fabsf((x1-x2)/(x1+x2));
    return fabsf(x1-x2);
  } else {
    return 1.0f;
  }
}

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

int main(void) {
    _f32 f, g1, g2, g3;
    unsigned int i, j, nerrs = 0;
    float diff, maxdiff = 0.0f;

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

    nerrs = 0;
    maxdiff = 0.0f;
    f.i = 0x3f800000;
    printf("start: %.7f\n", f.f-1.0f);
    for (i = 0; i < (1 << 23U); i++) {
        f.i = 0x3f800000|i;
        //f.f -= 1.0f;
        //g1.f = rcp(f.f);
        //g2.f = 1.0f/f.f;
        g1.f = fft_cbrtf(f.f);
        g2.f = 1.0f/cbrtf(f.f);
        if ((diff = rel_err_diff(g2.f, g1.f)) >=  0x1.0p-22f) {
            nerrs++;
            g3.f = cbrtf_estimate(f.f);
            if (diff > fabsf(maxdiff)) maxdiff = diff;
            if (!(i & 0xffff))
                printf("Error: rcp(%.7g): %.7g != 1.0f/(%.7g): %.7g (diff: %.7e / %a), g1: 0x%08x, g3: 0x%08x, g2: 0x%08x\n", f.f, g1.f, f.f, g2.f, diff, diff, g1.i, g3.i, g2.i);
        }
    }
    printf("Total errors: %u, largest error: cos: %.7e\n", nerrs, maxdiff);

    return 0;
}