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

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

static inline int fft_ilogbf(float x)
{
    _f32 hx;
    hx.f = fabsf(x);
    return ((hx.i >> 23) - 127);
}

static inline float ldexpf1(int n)
{
    _f32 tmp;
    /* normal result */
    n += 0x7f;
    tmp.i = (n << 23);
    return tmp.f;
}

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

float fft_cbrtf(float x)
{
    const float c0 = 0.333333333f;
    _f32 tmp;
    float r, t;
    uint32_t sign;

    tmp.f = x;
    sign = (tmp.i & 0x80000000);
    tmp.i &= 0x7fffffff;               /* x <- |x| */
    if(tmp.i==0||tmp.i>=0x7f800000)
        return(x);      /* cbrt(NaN,0,INF) is itself */

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

    /* new cbrt to 23 bits */
    //r  = (1.0f - x*t*t*t);
    //t += c0*t*(r - 2*c0*r*r);
    //r  = c0*(1.0f - x*t*t*t);
    //t += t*(r - 2*r*r);
    //r = (1/3)*(1.0 - x*r)
    //x = x + x*r*[1.0 - 2*r]
    //r = c0*(1.0f - x*t*t*t);
    r = (c0 - (c0*x)*t*t*t);
    t += t*(r + 2*r*r);
    r = (c0 - (c0*x)*t*t*t);
    t *= (1.0f + r + 2*r*r);
    //r = (1.0f - x*t*t*t);
    //t *= c0*(3.0f + r*(1.0f+2*c0*r));

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

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

#include <stdio.h>
#define ULPDIFF(given, expected) \
    (fabsf ((given) - (expected)) * ldexpf1(23 - fft_ilogbf(expected)))

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)
{
    _f32 f, g1, g2;
    unsigned int i, j, nerrs = 0;
    float diff, maxdiff = 0.0f;

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

    nerrs = 0;
    maxdiff = 0.0f;

    for (j = 127; j < 151; j++) {
        f.i = (j<<23);
        printf("start: %.7f (j: %u)\n", f.f, j);
        for (i = 0; i < (1 << 23U); i++) {
            f.i = (j<<23)|i;
            g1.f = cbrtf(f.f);
            g2.f = fft_cbrtf(f.f);
            //g1.f *= g1.f;
            if ((diff = fabsf(g1.f - g2.f)) >=  0x1.0p-21f) {
                long ulps = (long)ULPDIFF(g1.f, g2.f);
                nerrs++;
                if (diff > fabsf(maxdiff)) maxdiff = diff;
                if (ulps > 2) {
                    printf("Error: cbrtf(%.7g): %.7g != fft_cbrtf(%.7g): %.7g (diff: %.7e / %a, ULPs: %ld)\n", f.f, g1.f, f.f, g2.f, diff, diff, ulps);
                }
            }
        }
        printf("end: %.7f\n", f.f);
    }
    printf("Total errors: %u, largest error: %.7e\n", nerrs, maxdiff);

    return 0;
}