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

#define FLT_NAN         (*(float*)(&_flt_nan__))
#define FLT_INF         (*(float*)(&_flt_inf_positive__))
#define IsNan(_value)   ((*(uint32_t*)&_value)==_flt_nan__)
#define IsInf(_value)   (((*(uint32_t*)&_value)&0x7fffffff) == _flt_inf_positive)

uint32_t _flt_nan__ = 0x7fffffffL;
uint32_t _flt_inf_positive__ = 0x7f800000L;
uint32_t _flt_inf_negative__ = 0xff800000L;

#define LOG2    (0.6931471805599453094172321214581765680755f)
#define SQRT2    (1.4142135623737f)

/**
 * @brief calculate the exponential.
 * @param [in] x   input value of exp(x)
 * @return the result of exp(x)
 * @note Taylor deployment of exp function.
 */
float _expf(float x) {
    int i, k, negative = 0;
    float a, e, prev;

    //  make x as -ln(2) < x < ln(2)
    k = (int)(x / LOG2 + (x >= 0) ? 0.5f : -0.5f);
    x -= k * LOG2;

    if (x < 0) {
        negative = 1;
        x = -x;
    }

    e = 1 + x;
    a = x;
    i = 2;
    do {
        a *= x / (i++);
        prev = e;
        e += a;
    } while (a > FLT_MIN);
    if (negative != 0)
        e = 1 / e;
    return _ldexpf(e, k);
}

/**
 * @brief personal version of ldexp
 * @param [in] x value x
 * @param [in] exp
 * @result   x * 2^exp
 */
float    _ldexpf(float x, int exp) {
    float    twos_exponent;
    if (exp == 0)
        return    x;
    exp += 127;
    if (exp >= 255 || exp <= 0)
        return    FLT_NAN;
    exp <<= 23;
    twos_exponent = *(float*)(&exp);
    return    x * twos_exponent;
}

/**
 * @brief separate the  single real value as fraction and exponent
 */
float _frexpf(float x, int *exp) {
    uint32_t    val = *(uint32_t*)(&x);
    int    e = (val & 0x7f800000L) >> 23;
    uint32_t    frac = val & 0x7fffffL;
    if (e == 0) {
        *exp = 0;
        return    *(float*)(&val);
    }
    else if (e == 255) {
        *exp = 0;
        return    FLT_NAN;
    }
    else {
        e -= 127;
        *exp = e + 1;
        uint32_t    result = frac | (val & 0x80000000L);
        result |= ((-1) + 127) << 23;
        return    *(float*)&result;
    }
}

/**
 * @brief personal version of fabs
 * @param [in] x value x
 * @return absolute value of x
 */
float    _fabsf(float x) {
    uint32_t    value = *(uint32_t*)(&x);
    value &= 0x7fffffffL;    //    make sign bit to zero
    return *(float*)(&value);
}


/**
 * @brief calculate ln
 * @param [in] x   input value of ln(x)
 * @return the result of ln(x)
 * @note Euler transform of Taylor deployment of log function.
 */
float    _logf(float x) {
    int    ix, k;
    float    u, u2, s, prev;
    if (x <= 0)
        return    FLT_NAN;    //  error value
    //  make x over than -1 and under than 1.0 to prevent emerging
    _frexpf(x / SQRT2, &k);   //  k = (int)(log2(y)-0.5) (we decrease the exponent by 0.5 to make easy to converge.)
    x /= _ldexpf(1.0f, k);    //  x = x / 2^k;
    u = (x - 1) / (x + 1);    //  get log u with u = (x - 1) / (x + 1)
    u2 = u * u;
    ix = 1;
    s = u;
    do {
        u *= u2;
        ix += 2;
        prev = s;
        s += u / ix;
    } while (_fabsf(prev - s) > FLT_MIN);
    return LOG2 * k + 2 * s;
}


/**
* @brief  calculate power of x^y
* @param [in] x  number to be raised
* @param [in] y  power of x^y 
* @return   y th power of x
*/
float _powf(float x, float y) {

    if (y == (int)y) {
        float r = 1;
        int n = (int)((y > 0) ? y : (-y));
        while (n != 0) {
            if (n & 1) r *= x;
            x *= x;
            n >>= 1;
        }
        if (y < 0)
            r = 1 / r;
        return    r;
    }
    if (x > 0) {
        return    _expf(y*_logf(x));
    }
    else if (x != 0 || y <= 0) {
        return FLT_NAN;
    }
    return    0;
}

#if defined(_WIN32) || defined(_WIN64)
/**
* @brief  Test program
*/
int main() {
    float    val;
    float    frac;
    for (int j = 0; j < 5; ++j) {
        for (int i = 0; i <= 20; ++i) {
            val = (float)i / 3.3f;
            frac = (float)j / 2.2f;
            printf("%.2f : %.2f : %.4f : %.4f\n", frac, val, powf(frac, val), _powf(frac, val));
        }
    }
    getchar();

}
#endif