math.c

#include <stdint.h>

#define C_INF   (double)(1e+300 * 1e+300)
#define C_NAN   (double)(C_INF * 0.0)

// probably faster square root
// 0x1FF7A3BEA91D9B00   --  0.5*4503599627370496*(1023-0.0450465)
double C_sqrt(double num)
{
    // number should be positive
    if(num < 0.0) {
        // sqrt(< 0.0) can't be
        return C_NAN;
    }

    const double nhalf = (num * 0.5);   // needed below

    // See https://en.wikipedia.org/wiki/Fast_inverse_square_root
    int64_t i = *((int64_t *)&num);
    i = (0x1FF7A3BEA91D9B00 + (i >> 1));
    num = *((double *)&i);

    // lets calculate more precise value
    num *= (0.5 + (nhalf / (num * num)));
    num *= (0.5 + (nhalf / (num * num)));   // can be removed
    num *= (0.5 + (nhalf / (num * num)));   // can be removed
    return num;
}

// probably faster inverse square root
// 0x5FE6EC85E7DE30DA
// 0x5FE6EB3BFB58D000   --  1.5*4503599627370496*(1023-0.0450465)
double C_rsqrt(double num)
{
    // number should be positive non-zero
    if(num <= 0.0) {
        return (num ? C_NAN : C_INF);
    }

    const double nhalf = (num * 0.5);   // needed below

    // See https://en.wikipedia.org/wiki/Fast_inverse_square_root
    int64_t i = *((int64_t *)&num);
    i = (0x5FE6EB3BFB58D000 - (i >> 1));
    num = *((double *)&i);

    // lets calculate more precise value
    num *= (1.5 - (nhalf * num * num));
    num *= (1.5 - (nhalf * num * num)); // can be removed
    num *= (1.5 - (nhalf * num * num)); // can be removed
    return num;
}