#include <stdint.h>
#include <string.h>
#include <math.h>
//0.28BE60DB9391054A7F09D5F47D4D377036D8A5664F10E4107F9458EB
//0.517CC1B727220A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D6

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

/* 190 bits of 2/pi for Payne-Hanek style argument reduction. */
static const unsigned int two_over_pi_f [] = 
{
    0x00000000,
    0x28be60db,
    0x9391054a,
    0x7f09d5f4,
    0x7d4d3770,
    0x36d8a566,
    0x4f10e411
};

float mod_2pi_slow (float a)
{
    unsigned long int p, ia;
    unsigned int hi, mid, lo, tmp, i;
    int e;
    double d;
    _f32 tmp2;
    float r;

    a = fabsf(a);
    tmp2.f = a;
    //if (a <= 0x1.921fb5p+02f || a != a) {
    if (tmp2.i <= 0x40490fda || a != a) {
        return a;
    }
    e = (int)(tmp2.i>>23)-127;
    tmp2.i &= 0x007fffff;
    //tmp2.i |= 0x4f000000;
    //ia = (unsigned int)(tmp2.f);
    ia = (unsigned int)(tmp2.i<<8)|0x80000000;

    /* extract 96 relevant bits of 2/pi based on magnitude of argument */ 
    i = (e < 0) ? 0 : (unsigned int)e >> 5;
    e = (unsigned int)e & 31;

    hi  = two_over_pi_f [i+0];
    mid = two_over_pi_f [i+1];
    lo  = two_over_pi_f [i+2];
    tmp = two_over_pi_f [i+3];

    if (e) {
        hi  = (hi  << e) | (mid >> (32 - e));
        mid = (mid << e) | (lo  >> (32 - e));
        lo  = (lo  << e) | (tmp >> (32 - e));
    }

    /* compute product x * 2/pi in 2.62 fixed-point format */
    p = ia * lo;
    p = ia * mid + (p >> 32);
    p = ((ia * hi) << 32) + p;

    /* round quotient to nearest */
    p <<= 1;
    p >>= 1;
    //p = p & 0x7fffffffffffffffULL;     // fraction

    /* compute remainder of x * (2/pi) */
    //d = (double)p;
    //d = d * 0x1.921fb54442d18p+0; // 1.5707963267948966
    //r = (float)d * 0x1.0p-61f;
    r = (float)p * 0x1.0p-61f;
    r = r * 0x1.921fb54442d18p+00f; // 1.5707963267948966 * 0x1.0p-61
    return r;
}