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

#include "core_math_private.h"

void QuatMultiply(const vec4_t left, const vec4_t right, vec4_t out)
{
#if EMULATE_X87
    /* x87-faithful transcription (stock coduo_lnxded 0x0806960c): each component
     * is a whole-expression 80-bit sum/difference of four lane products, rounded
     * to float only at the store — stock keeps the intermediates in x87 registers,
     * so the multi-op chains run in the shim. */
    out[0] = x87f_store_f32(x87f_sub(
        x87f_add(x87f_add(x87f_mul(x87f_load_f32(left[0]), x87f_load_f32(right[3])),
                          x87f_mul(x87f_load_f32(left[3]), x87f_load_f32(right[0]))),
                 x87f_mul(x87f_load_f32(left[2]), x87f_load_f32(right[1]))),
        x87f_mul(x87f_load_f32(left[1]), x87f_load_f32(right[2]))));
    out[1] = x87f_store_f32(x87f_add(
        x87f_add(x87f_sub(x87f_mul(x87f_load_f32(left[1]), x87f_load_f32(right[3])),
                          x87f_mul(x87f_load_f32(left[2]), x87f_load_f32(right[0]))),
                 x87f_mul(x87f_load_f32(left[3]), x87f_load_f32(right[1]))),
        x87f_mul(x87f_load_f32(left[0]), x87f_load_f32(right[2]))));
    out[2] = x87f_store_f32(x87f_add(
        x87f_sub(x87f_add(x87f_mul(x87f_load_f32(left[2]), x87f_load_f32(right[3])),
                          x87f_mul(x87f_load_f32(left[1]), x87f_load_f32(right[0]))),
                 x87f_mul(x87f_load_f32(left[0]), x87f_load_f32(right[1]))),
        x87f_mul(x87f_load_f32(left[3]), x87f_load_f32(right[2]))));
    out[3] = x87f_store_f32(x87f_sub(
        x87f_sub(x87f_sub(x87f_mul(x87f_load_f32(left[3]), x87f_load_f32(right[3])),
                          x87f_mul(x87f_load_f32(left[0]), x87f_load_f32(right[0]))),
                 x87f_mul(x87f_load_f32(left[1]), x87f_load_f32(right[1]))),
        x87f_mul(x87f_load_f32(left[2]), x87f_load_f32(right[2]))));
#else
    out[0] = (((left[0] * right[3]) + (left[3] * right[0])) +
              (left[2] * right[1])) -
             (left[1] * right[2]);
    out[1] = (((left[1] * right[3]) - (left[2] * right[0])) +
              (left[3] * right[1])) +
             (left[0] * right[2]);
    out[2] = (((left[2] * right[3]) + (left[1] * right[0])) -
              (left[0] * right[1])) +
             (left[3] * right[2]);
    out[3] = (((left[3] * right[3]) - (left[0] * right[0])) -
              (left[1] * right[1])) -
             (left[2] * right[2]);
#endif
}
