// ---------------------------- VSISQRT.C -----------------------
// VSpace 3d library.
// Designed and written by Javier Ar‚valo Baeza.
// Adaptation of the Fast inverse square root algorithm from Ken
// Turkowski, Graphics Gems V
// And the square root algo from Gems 1 by Paul Lalonde and Robert Dawson
// Make sure the float parameters are really float and not double.

#ifdef SQRTAPP

#include "vsisqrt.h"

#include <math.h>

    // Define lookup table.
long VSISQ_LookupTable[VSISQ_TABLE_SIZE];
long VSISQ_SqrtTab[0x100];


typedef union {
    dword i;
    float f;
} _flint;

static _flint fi, fo;

static void build_table(void) {
    unsigned short i;
    float f;
    unsigned int *fi = (unsigned*)&f;   // To access the bits of a float in
                                        // C quickly we must misuse pointers
    for (i = 0; i <= 0x7f; i++) {
        *fi = 0;

            // Build a float with the bit pattern i as mantissa
            // and an exponent of 0, stored as 127
        *fi = (i << 16) | (127 << 23);
        f = sqrt(f);

            // Take the square root then strip the first 7 bits of
            // the mantissa into the table
        VSISQ_SqrtTab[i] = (*fi & 0x7fffff);

            // Repeat the process, this time with an exponent of 1,
            // stored as 128
        *fi = 0;
        *fi = (i << 16) | (128 << 23);
        f = sqrt(f);
        VSISQ_SqrtTab[i+0x80] = (*fi & 0x7fffff);
    }
}


extern void VSISQ_Init(void) {
    long f;
    long *h;

    for (f = 0, h = VSISQ_LookupTable; f < VSISQ_TABLE_SIZE; f++) {
        fi.i = ((VSISQ_EXP_BIAS-1) << VSISQ_EXP_POS) | (f << VSISQ_LOOKUP_POS);
        fo.f = 1.0 / sqrt(fi.f);
        *h++ = ((fo.i + (1 << (VSISQ_SEED_POS-2))) >> VSISQ_SEED_POS) & 0x007FFFFF;
    }
    VSISQ_LookupTable[VSISQ_TABLE_SIZE/2] = 0x007FFFFF;   // Special case for 1.0
    build_table();
}

/*
extern VSISQ_TFloat VSISQ_ISqrt(float x) {
    dword a = ((_flint*)(&x))->i;
    float arg = x;
    static _flint seed;
    VSISQ_TFloat r;


//    seed.i = VSISQ_SET_EXP(((3*VSISQ_EXP_BIAS-1) - VSISQ_GET_EXP(a)) >> 1)
//           | VSISQ_SET_MANTSEED(VSISQ_LookupTable[VSISQ_GET_EMANT(a)]);

//    In other words:
    seed.i = (((((3*127-1) << 22)) - ((a & 0x7F800000) >> 1)) & 0x7F800000)
           | VSISQ_LookupTable[(a >> 15) & 0x1FF];

        // Init seed, accurate to LOOKUP_BITS
    r = seed.f;
        // First iteration, accurate to 2*LOOKUP_BITS
    r = (3.0 - r * r * arg) * r * 0.5;
        // Second iteration, accurate to 4*LOOKUP_BITS
//    r = (3.0 - r * r * arg) * r * 0.5;
        // Third iteration, accurate to 8*LOOKUP_BITS
//    r = (3.0 - r * r * arg) * r * 0.5;
    return r;
}
*/

#endif
// ---------------------------- VSISQRT.C -----------------------
