#include <math.h>

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

#include "core_math_private.h"

long double Distance(const vec3_t left, const vec3_t right)
{
    float length;

#if EMULATE_X87
    /* x87-faithful transcription (stock coduo_lnxded 0x0806657b): deltas
     * right[i]-left[i] rounded to float slots (fstps), squares summed
     * left-to-right in 80-bit (faddp), the sum rounded straight to double64 for
     * the CRT sqrt (fstpl — no float spill), sqrt narrowed to float and returned
     * via the long-double ABI. sqrt is IEEE correctly-rounded (bit-portable). */
    const float x =
        x87f_store_f32(x87f_sub(x87f_load_f32(right[0]), x87f_load_f32(left[0])));
    const float y =
        x87f_store_f32(x87f_sub(x87f_load_f32(right[1]), x87f_load_f32(left[1])));
    const float z =
        x87f_store_f32(x87f_sub(x87f_load_f32(right[2]), x87f_load_f32(left[2])));
    x87f sum = x87f_mul(x87f_load_f32(x), x87f_load_f32(x));
    sum = x87f_add(sum, x87f_mul(x87f_load_f32(y), x87f_load_f32(y)));
    sum = x87f_add(sum, x87f_mul(x87f_load_f32(z), x87f_load_f32(z)));
    length = (float)sqrt(x87f_store_f64(sum));
#else
    const float x = right[0] - left[0];
    const float y = right[1] - left[1];
    const float z = right[2] - left[2];
    length = (float)sqrt((double)(x * x + y * y + z * z));
#endif
    return length;
}
