#include "coduo_x87emu.h" /* defines EMULATE_X87; x87 shim when it is 1 */

#include "core_math_private.h"

long double QuatNormalizedVectorPartLengthSquared(const vec4_t vector)
{
    float xSquared;
    float ySquared;
    float zSquared;
    float denominator;
    float result;

#if EMULATE_X87
    /* x87-faithful transcription (stock coduo_lnxded 0x0806998f): each squared
     * component is rounded to its float slot; the denominator sum
     * ((x^2+y^2)+z^2)+(w*w) is accumulated in 80-bit and rounded to float; when
     * nonzero the reciprocal, the three scaled squares, and their summed result
     * are each computed in 80-bit and rounded to float. */
    xSquared = x87f_store_f32(x87f_mul(x87f_load_f32(vector[0]),
                                       x87f_load_f32(vector[0])));
    ySquared = x87f_store_f32(x87f_mul(x87f_load_f32(vector[1]),
                                       x87f_load_f32(vector[1])));
    zSquared = x87f_store_f32(x87f_mul(x87f_load_f32(vector[2]),
                                       x87f_load_f32(vector[2])));
    denominator = x87f_store_f32(x87f_add(
        x87f_add(x87f_add(x87f_load_f32(xSquared), x87f_load_f32(ySquared)),
                 x87f_load_f32(zSquared)),
        x87f_mul(x87f_load_f32(vector[3]), x87f_load_f32(vector[3]))));

    if (denominator != 0.0f) {
        denominator = x87f_store_f32(
            x87f_div(x87f_load_f32(1.0f), x87f_load_f32(denominator)));
        xSquared = x87f_store_f32(
            x87f_mul(x87f_load_f32(xSquared), x87f_load_f32(denominator)));
        ySquared = x87f_store_f32(
            x87f_mul(x87f_load_f32(ySquared), x87f_load_f32(denominator)));
        zSquared = x87f_store_f32(
            x87f_mul(x87f_load_f32(zSquared), x87f_load_f32(denominator)));
        result = x87f_store_f32(
            x87f_add(x87f_add(x87f_load_f32(xSquared), x87f_load_f32(ySquared)),
                     x87f_load_f32(zSquared)));
    } else {
        result = 0.0f;
    }
#else
    xSquared = vector[0] * vector[0];
    ySquared = vector[1] * vector[1];
    zSquared = vector[2] * vector[2];
    denominator = ((xSquared + ySquared) + zSquared) + (vector[3] * vector[3]);

    if (denominator != 0.0f) {
        denominator = 1.0f / denominator;
        xSquared = xSquared * denominator;
        ySquared = ySquared * denominator;
        zSquared = zSquared * denominator;
        result = (xSquared + ySquared) + zSquared;
    } else {
        result = 0.0f;
    }
#endif

    return (long double)result;
}
