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

#include "core_math_private.h"
#include "coduo_ctype_compat.h"
#include "coduo_int32_bits.h"
#include "coduo_native_x87.h"

#include <ctype.h>
#include <math.h>
#include <stddef.h>
#include <string.h>

enum {
    CODUO_COLOR_BYTE_SCALE = 255,
    CODUO_COLOR_ALPHA_OPAQUE = 0xff,
    CODUO_ANGLE_PACK_SCALE = 65536,
    CODUO_ANGLE_PACK_MASK = 0xffff,
    CODUO_WILDCARD_CHUNK_SIZE = 1024,
    CODUO_NORMAL_ENCODE_SCALE = 255,
    CODUO_NORMAL_BYTE_SCALE = 360,
    CODUO_NORMAL_BYTE_RANGE = 256,
    CODUO_NORMAL_BYTE_NEGATIVE_Z = 0x80,
    CODUO_PATH_MATCH_BUFFER_SIZE = 64,
    CODUO_HASH_POSITION_BIAS = 119
};

#define CODUO_ANGLE_180 180.0f
#define CODUO_ANGLE_360 360.0f
#define CODUO_MILLISECONDS_PER_SECOND 1000.0f

uint32_t randSeed = 0x89abcdefU;

uint32_t ColorBytes3(float red, float green, float blue)
{
    uint8_t bytes[4];

#if EMULATE_X87
    /* red*scale is an 80-bit product truncated toward zero to int (stock keeps it
     * in a register before fistp, no float spill). */
    bytes[0] = (uint8_t)(int16_t)x87f_store_i32_trunc(x87f_mul(
        x87f_load_f32(red), x87f_load_f32((float)CODUO_COLOR_BYTE_SCALE)));
    bytes[1] = (uint8_t)(int16_t)x87f_store_i32_trunc(x87f_mul(
        x87f_load_f32(green), x87f_load_f32((float)CODUO_COLOR_BYTE_SCALE)));
    bytes[2] = (uint8_t)(int16_t)x87f_store_i32_trunc(x87f_mul(
        x87f_load_f32(blue), x87f_load_f32((float)CODUO_COLOR_BYTE_SCALE)));
#elif defined(__x86_64__)
    static const float scale = (float)CODUO_COLOR_BYTE_SCALE;
    coduo_x87_truncation_control_t control;
    int16_t converted[3];

    converted[0] = CODUO_X87_TRUNCATE_I16_FIRST(
        &control, (long double)red * (long double)scale);
    converted[1] = CODUO_X87_TRUNCATE_I16_NEXT(
        &control, (long double)green * (long double)scale);
    converted[2] = CODUO_X87_TRUNCATE_I16_NEXT(
        &control, (long double)blue * (long double)scale);

    bytes[0] = (uint8_t)converted[0];
    bytes[1] = (uint8_t)converted[1];
    bytes[2] = (uint8_t)converted[2];
#else
    bytes[0] = (uint8_t)(int16_t)(red * (float)CODUO_COLOR_BYTE_SCALE);
    bytes[1] = (uint8_t)(int16_t)(green * (float)CODUO_COLOR_BYTE_SCALE);
    bytes[2] = (uint8_t)(int16_t)(blue * (float)CODUO_COLOR_BYTE_SCALE);
#endif
    bytes[3] = CODUO_COLOR_ALPHA_OPAQUE;

    uint32_t packed;
    memcpy(&packed, bytes, sizeof(packed));
    return packed;
}

uint32_t ColorBytes4(float red, float green, float blue, float alpha)
{
    uint8_t bytes[4];

#if EMULATE_X87
    bytes[0] = (uint8_t)(int16_t)x87f_store_i32_trunc(x87f_mul(
        x87f_load_f32(red), x87f_load_f32((float)CODUO_COLOR_BYTE_SCALE)));
    bytes[1] = (uint8_t)(int16_t)x87f_store_i32_trunc(x87f_mul(
        x87f_load_f32(green), x87f_load_f32((float)CODUO_COLOR_BYTE_SCALE)));
    bytes[2] = (uint8_t)(int16_t)x87f_store_i32_trunc(x87f_mul(
        x87f_load_f32(blue), x87f_load_f32((float)CODUO_COLOR_BYTE_SCALE)));
    bytes[3] = (uint8_t)(int16_t)x87f_store_i32_trunc(x87f_mul(
        x87f_load_f32(alpha), x87f_load_f32((float)CODUO_COLOR_BYTE_SCALE)));
#elif defined(__x86_64__)
    static const float scale = (float)CODUO_COLOR_BYTE_SCALE;
    coduo_x87_truncation_control_t control;
    int16_t converted[4];

    converted[0] = CODUO_X87_TRUNCATE_I16_FIRST(
        &control, (long double)red * (long double)scale);
    converted[1] = CODUO_X87_TRUNCATE_I16_NEXT(
        &control, (long double)green * (long double)scale);
    converted[2] = CODUO_X87_TRUNCATE_I16_NEXT(
        &control, (long double)blue * (long double)scale);
    converted[3] = CODUO_X87_TRUNCATE_I16_NEXT(
        &control, (long double)alpha * (long double)scale);

    bytes[0] = (uint8_t)converted[0];
    bytes[1] = (uint8_t)converted[1];
    bytes[2] = (uint8_t)converted[2];
    bytes[3] = (uint8_t)converted[3];
#else
    bytes[0] = (uint8_t)(int16_t)(red * (float)CODUO_COLOR_BYTE_SCALE);
    bytes[1] = (uint8_t)(int16_t)(green * (float)CODUO_COLOR_BYTE_SCALE);
    bytes[2] = (uint8_t)(int16_t)(blue * (float)CODUO_COLOR_BYTE_SCALE);
    bytes[3] = (uint8_t)(int16_t)(alpha * (float)CODUO_COLOR_BYTE_SCALE);
#endif

    uint32_t packed;
    memcpy(&packed, bytes, sizeof(packed));
    return packed;
}

long double NormalizeColor(const vec3_t source, vec3_t dest)
{
    float max = source[0];

    if (max < source[1]) {
        max = source[1];
    }
    if (max < source[2]) {
        max = source[2];
    }

    if (max == 0.0f) {
        dest[0] = 0.0f;
        dest[1] = 0.0f;
        dest[2] = 0.0f;
    } else {
#if EMULATE_X87
        /* Each source[i]/max is a single 80-bit x87 divide rounded to float
         * (fdivrs/fstps); the max-selection compares are plain float compares. */
        dest[0] = x87f_store_f32(
            x87f_div(x87f_load_f32(source[0]), x87f_load_f32(max)));
        dest[1] = x87f_store_f32(
            x87f_div(x87f_load_f32(source[1]), x87f_load_f32(max)));
        dest[2] = x87f_store_f32(
            x87f_div(x87f_load_f32(source[2]), x87f_load_f32(max)));
#else
        dest[0] = source[0] / max;
        dest[1] = source[1] / max;
        dest[2] = source[2] / max;
#endif
    }

    return max;
}

long double AngleMod(float angleDegrees)
{
#if EMULATE_X87
    /* x87-faithful transcription (stock coduo_lnxded 0x08069da8): angleDegrees *
     * (65536/360) is formed in 80-bit and truncated toward zero to int (fistp
     * under a truncate control word); the masked integer is reloaded (fild) and
     * multiplied by (360/65536) in 80-bit. Stock returns that product unrounded
     * in st0; no engine caller consumes it (decl only), so it is rounded to float
     * here to match the float-projected digest and avoid an arm64 truncation. */
    int32_t packed = x87f_store_i32_trunc(
        x87f_mul(x87f_load_f32(angleDegrees),
                 x87f_load_f32((float)CODUO_ANGLE_PACK_SCALE / CODUO_ANGLE_360)));
    return (long double)x87f_store_f32(x87f_mul(
        x87f_load_f32(CODUO_ANGLE_360 / (float)CODUO_ANGLE_PACK_SCALE),
        x87f_load_i32(packed & CODUO_ANGLE_PACK_MASK)));
#elif defined(__x86_64__)
    static const float packScale =
        (float)CODUO_ANGLE_PACK_SCALE / CODUO_ANGLE_360;
    int32_t packed = CODUO_X87_TRUNCATE_I32(
        (long double)angleDegrees * (long double)packScale);

    return (float)(CODUO_ANGLE_360 / (float)CODUO_ANGLE_PACK_SCALE) *
           (float)(packed & CODUO_ANGLE_PACK_MASK);
#else
    int32_t packed =
        (int32_t)(angleDegrees * ((float)CODUO_ANGLE_PACK_SCALE /
                                  CODUO_ANGLE_360));

    return (float)(CODUO_ANGLE_360 / (float)CODUO_ANGLE_PACK_SCALE) *
           (float)(packed & CODUO_ANGLE_PACK_MASK);
#endif
}

long double LerpAngle(float from, float to, float frac)
{
#if EMULATE_X87
    /* x87-faithful transcription (stock coduo_lnxded 0x08069deb): each `to - from`
     * (including the two branch tests) is the RAW 80-bit difference — stock never
     * spills it to a float before the fucompp, so the comparisons must run in
     * 80-bit; the 360 wraps and the final from + (to-from)*frac are 80-bit and
     * rounded to float. */
    {
        x87f diff = x87f_sub(x87f_load_f32(to), x87f_load_f32(from));
        if (x87f_lt(x87f_load_f32(CODUO_ANGLE_180), diff)) { /* diff > 180 */
            to = x87f_store_f32(
                x87f_sub(x87f_load_f32(to), x87f_load_f32(CODUO_ANGLE_360)));
        }
        diff = x87f_sub(x87f_load_f32(to), x87f_load_f32(from));
        if (x87f_lt(diff, x87f_load_f32(-CODUO_ANGLE_180))) { /* diff < -180 */
            to = x87f_store_f32(
                x87f_add(x87f_load_f32(to), x87f_load_f32(CODUO_ANGLE_360)));
        }
        {
            x87f d = x87f_sub(x87f_load_f32(to), x87f_load_f32(from));
            float result = x87f_store_f32(
                x87f_add(x87f_mul(d, x87f_load_f32(frac)), x87f_load_f32(from)));
            return result;
        }
    }
#else
    if (to - from > CODUO_ANGLE_180) {
        to -= CODUO_ANGLE_360;
    }
    if (to - from < -CODUO_ANGLE_180) {
        to += CODUO_ANGLE_360;
    }

    /* stock 0x8069e4d rounds the result to float before the x87 return */
    float result = from + ((to - from) * frac);
    return result;
#endif
}

long double AngleSubtract(float angle1, float angle2)
{
    /* ORIGINAL_BINARY_BUG: infinity or a sufficiently large finite binary32
     * value can make a 360-degree store leave delta unchanged, so either
     * loop can fail to terminate; see
     * original-bugs/client-loop-angle-normalizers-nontermination.md. */
#if EMULATE_X87
    /* delta and each 360 wrap are 80-bit fsub/fadd rounded to float; loop
     * compares are plain float compares (stock coduo_lnxded 0x08069e5b). */
    float delta =
        x87f_store_f32(x87f_sub(x87f_load_f32(angle1), x87f_load_f32(angle2)));

    while (delta > CODUO_ANGLE_180) {
        delta = x87f_store_f32(
            x87f_sub(x87f_load_f32(delta), x87f_load_f32(CODUO_ANGLE_360)));
    }
    while (delta < -CODUO_ANGLE_180) {
        delta = x87f_store_f32(
            x87f_add(x87f_load_f32(delta), x87f_load_f32(CODUO_ANGLE_360)));
    }
#else
    float delta = angle1 - angle2;

    while (delta > CODUO_ANGLE_180) {
        delta -= CODUO_ANGLE_360;
    }
    while (delta < -CODUO_ANGLE_180) {
        delta += CODUO_ANGLE_360;
    }
#endif

    return delta;
}

void AnglesSubtract(const vec3_t angles1, const vec3_t angles2, vec3_t out)
{
    out[0] = (float)AngleSubtract(angles1[0], angles2[0]);
    out[1] = (float)AngleSubtract(angles1[1], angles2[1]);
    out[2] = (float)AngleSubtract(angles1[2], angles2[2]);
}

long double AngleNormalize360(float angleDegrees)
{
#if EMULATE_X87
    /* Identical packing to AngleMod (stock coduo_lnxded 0x08069f27): 80-bit
     * angleDegrees*(65536/360) truncated to int, masked, reloaded and multiplied
     * by (360/65536) in 80-bit. Stock returns that product unrounded in st0;
     * this existing non-x87 path retains its permitted binary32 projection. */
    int32_t packed = x87f_store_i32_trunc(
        x87f_mul(x87f_load_f32(angleDegrees),
                 x87f_load_f32((float)CODUO_ANGLE_PACK_SCALE / CODUO_ANGLE_360)));
    return (long double)x87f_store_f32(x87f_mul(
        x87f_load_f32(CODUO_ANGLE_360 / (float)CODUO_ANGLE_PACK_SCALE),
        x87f_load_i32(packed & CODUO_ANGLE_PACK_MASK)));
#elif defined(__x86_64__)
    static const float packScale =
        (float)CODUO_ANGLE_PACK_SCALE / CODUO_ANGLE_360;
    int32_t packed = CODUO_X87_TRUNCATE_I32(
        (long double)angleDegrees * (long double)packScale);

    return (float)(CODUO_ANGLE_360 / (float)CODUO_ANGLE_PACK_SCALE) *
           (float)(packed & CODUO_ANGLE_PACK_MASK);
#else
    int32_t packed =
        (int32_t)(angleDegrees * ((float)CODUO_ANGLE_PACK_SCALE /
                                  CODUO_ANGLE_360));

    return (float)(CODUO_ANGLE_360 / (float)CODUO_ANGLE_PACK_SCALE) *
           (float)(packed & CODUO_ANGLE_PACK_MASK);
#endif
}

long double AngleNormalize180(float angleDegrees)
{
    float angle = (float)AngleNormalize360(angleDegrees);

#if EMULATE_X87
    /* the single 360 subtract is an 80-bit fsub rounded to float; the compare is
     * a plain float compare (stock coduo_lnxded 0x08069f6a). */
    if (angle > CODUO_ANGLE_180) {
        angle = x87f_store_f32(
            x87f_sub(x87f_load_f32(angle), x87f_load_f32(CODUO_ANGLE_360)));
    }
#else
    if (angle > CODUO_ANGLE_180) {
        angle -= CODUO_ANGLE_360;
    }
#endif

    return angle;
}

long double AngleNormalize360Accurate(float angleDegrees)
{
    /* ORIGINAL_BINARY_BUG: each adjusted value is stored to binary32.
     * Infinity or a sufficiently large finite value can therefore make the
     * selected loop fail to progress; see
     * original-bugs/client-loop-angle-normalizers-nontermination.md.
     *
     * ORIGINAL_BINARY_BUG: the negative arm returns without entering the
     * subtract arm. A tiny negative value whose addition rounds to 360.0f is
     * therefore returned as 360.0f; see
     * original-bugs/client-angle-normalize360-negative-boundary.md. */
    if (angleDegrees < 0.0f) {
#if EMULATE_X87
        /* each wrap is an 80-bit operation rounded to float; comparisons are
         * plain float comparisons (stock coduo_lnxded 0x08069fab). */
        do {
            angleDegrees = x87f_store_f32(
                x87f_add(x87f_load_f32(angleDegrees),
                         x87f_load_f32(CODUO_ANGLE_360)));
        } while (angleDegrees < 0.0f);
#else
        do {
            angleDegrees += CODUO_ANGLE_360;
        } while (angleDegrees < 0.0f);
#endif
    } else if (angleDegrees >= CODUO_ANGLE_360) {
#if EMULATE_X87
        do {
            angleDegrees = x87f_store_f32(
                x87f_sub(x87f_load_f32(angleDegrees),
                         x87f_load_f32(CODUO_ANGLE_360)));
        } while (angleDegrees >= CODUO_ANGLE_360);
#else
        do {
            angleDegrees -= CODUO_ANGLE_360;
        } while (angleDegrees >= CODUO_ANGLE_360);
#endif
    }

    return angleDegrees;
}

long double AngleNormalize180Accurate(float angleDegrees)
{
    /* ORIGINAL_BINARY_BUG: infinity or a sufficiently large finite binary32
     * value can make a 360-degree store leave angleDegrees unchanged, so
     * either loop can fail to terminate; see
     * original-bugs/client-loop-angle-normalizers-nontermination.md. */
#if EMULATE_X87
    /* each 360 wrap is an 80-bit fadd/fsub rounded to float (stock coduo_lnxded
     * 0x0806a028). */
    if (angleDegrees <= -CODUO_ANGLE_180) {
        do {
            angleDegrees = x87f_store_f32(x87f_add(
                x87f_load_f32(angleDegrees), x87f_load_f32(CODUO_ANGLE_360)));
        } while (angleDegrees <= -CODUO_ANGLE_180);
    } else if (angleDegrees > CODUO_ANGLE_180) {
        do {
            angleDegrees = x87f_store_f32(x87f_sub(
                x87f_load_f32(angleDegrees), x87f_load_f32(CODUO_ANGLE_360)));
        } while (angleDegrees > CODUO_ANGLE_180);
    }
#else
    if (angleDegrees <= -CODUO_ANGLE_180) {
        do {
            angleDegrees += CODUO_ANGLE_360;
        } while (angleDegrees <= -CODUO_ANGLE_180);
    } else if (angleDegrees > CODUO_ANGLE_180) {
        do {
            angleDegrees -= CODUO_ANGLE_360;
        } while (angleDegrees > CODUO_ANGLE_180);
    }
#endif

    return angleDegrees;
}

long double AngleDelta(float angle1, float angle2)
{
#if EMULATE_X87
    /* angle1-angle2 is an 80-bit fsub rounded to the float argument passed to
     * AngleNormalize180 (stock coduo_lnxded 0x0806a0ad). */
    return AngleNormalize180(
        x87f_store_f32(x87f_sub(x87f_load_f32(angle1), x87f_load_f32(angle2))));
#else
    return AngleNormalize180(angle1 - angle2);
#endif
}

float RadiusFromBounds(const vec3_t mins, const vec3_t maxs)
{
    vec3_t corner;

    for (int32_t axis = 0; axis < 3; ++axis) {
        float min = fabsf(mins[axis]);
        float max = fabsf(maxs[axis]);
        corner[axis] = min > max ? min : max;
    }

#if EMULATE_X87
    /* fabsf and the min>max selection are exact/plain-compare native ops; only the
     * squared sum through double64 sqrt to float is x87 (stock coduo_lnxded
     * 0x0806a0c3). */
    return (float)sqrt(x87f_store_f64(x87f_add(
        x87f_add(x87f_mul(x87f_load_f32(corner[0]), x87f_load_f32(corner[0])),
                 x87f_mul(x87f_load_f32(corner[1]), x87f_load_f32(corner[1]))),
        x87f_mul(x87f_load_f32(corner[2]), x87f_load_f32(corner[2])))));
#else
    return (float)sqrt(
        (double)((corner[0] * corner[0]) + (corner[1] * corner[1]) +
                 (corner[2] * corner[2])));
#endif
}

void YawToAxis(float angleDegrees, coduo_compact_affine_matrix43_t *matrix)
{
    vec3_t right;

    YawVectors(angleDegrees, matrix->axis[0], right);
    matrix->axis[2][0] = 0.0f;
    matrix->axis[2][1] = 0.0f;
    matrix->axis[2][2] = 1.0f;
    /* Stock 0x806a515..0x806a540 loads each vec3_origin lane and subtracts
     * right[i]. This is arithmetic positive-zero minus right, not a sign-bit
     * toggle, and it preserves the original global dependency and zero sign. */
    matrix->axis[1][0] = vec3_origin[0] - right[0];
    matrix->axis[1][1] = vec3_origin[1] - right[1];
    matrix->axis[1][2] = vec3_origin[2] - right[2];
}

void AxisToAngles(const vec3_t axis[3], vec3_t angles)
{
    vec3_t axis1;
    float sinValue;
    float cosValue;
    float yawTemp;

    vectoangles(axis[0], angles);

    memcpy(&axis1[0], &axis[1][0], sizeof(axis1[0]));
    memcpy(&axis1[1], &axis[1][1], sizeof(axis1[1]));
    memcpy(&axis1[2], &axis[1][2], sizeof(axis1[2]));

    /* stock multiplies by double pi then divides by double 180.0 (e.g.
     * 0x806a58d..0x806a59b); a folded pi/180 constant rounds differently */
#if EMULATE_X87
    SinCosFloat(
        x87f_store_f32(x87f_div(
            x87f_mul(x87f_neg(x87f_load_f32(angles[1])), x87f_load_f64(M_PI)),
            x87f_load_f64((double)CODUO_ANGLE_180))),
        &sinValue, &cosValue);
    yawTemp = x87f_store_f32(
        x87f_sub(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(axis1[0])),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(axis1[1]))));
    axis1[1] = x87f_store_f32(
        x87f_add(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(axis1[1])),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(axis1[0]))));

    SinCosFloat(
        x87f_store_f32(x87f_div(
            x87f_mul(x87f_neg(x87f_load_f32(angles[0])), x87f_load_f64(M_PI)),
            x87f_load_f64((double)CODUO_ANGLE_180))),
        &sinValue, &cosValue);
    axis1[0] = x87f_store_f32(
        x87f_add(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(yawTemp)),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(axis1[2]))));
    axis1[2] = x87f_store_f32(
        x87f_sub(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(axis1[2])),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(yawTemp))));
#else
    SinCosFloat((float)(((double)-angles[1] * M_PI) /
                           (double)CODUO_ANGLE_180),
                 &sinValue, &cosValue);
    yawTemp = (cosValue * axis1[0]) - (sinValue * axis1[1]);
    axis1[1] = (cosValue * axis1[1]) + (sinValue * axis1[0]);

    SinCosFloat((float)(((double)-angles[0] * M_PI) /
                           (double)CODUO_ANGLE_180),
                 &sinValue, &cosValue);
    axis1[0] = (cosValue * yawTemp) + (sinValue * axis1[2]);
    axis1[2] = (cosValue * axis1[2]) - (sinValue * yawTemp);
#endif

    float roll = (float)vectosignedpitch(axis1);
    if (axis1[1] < 0.0f) {
        if (roll < 0.0f) {
            roll += CODUO_ANGLE_180;
        } else {
            roll -= CODUO_ANGLE_180;
        }
        angles[2] = roll;
    } else {
        uint32_t rollBits;
        memcpy(&rollBits, &roll, sizeof(rollBits));
        rollBits ^= CODUO_FLOAT_SIGN_BIT;
        memcpy(&angles[2], &rollBits, sizeof(angles[2]));
    }
}

void Axis4ToAngles(const DObjSkelMat *matrix, vec3_t angles)
{
    vec3_t axis1;
    float sinValue;
    float cosValue;
    float yawTemp;

    vectoangles(matrix->axis[0], angles);

    memcpy(&axis1[0], &matrix->axis[1][0], sizeof(axis1[0]));
    memcpy(&axis1[1], &matrix->axis[1][1], sizeof(axis1[1]));
    memcpy(&axis1[2], &matrix->axis[1][2], sizeof(axis1[2]));

    /* stock multiplies by double pi then divides by double 180.0 (e.g.
     * 0x806a58d..0x806a59b); a folded pi/180 constant rounds differently */
#if EMULATE_X87
    SinCosFloat(
        x87f_store_f32(x87f_div(
            x87f_mul(x87f_neg(x87f_load_f32(angles[1])), x87f_load_f64(M_PI)),
            x87f_load_f64((double)CODUO_ANGLE_180))),
        &sinValue, &cosValue);
    yawTemp = x87f_store_f32(
        x87f_sub(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(axis1[0])),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(axis1[1]))));
    axis1[1] = x87f_store_f32(
        x87f_add(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(axis1[1])),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(axis1[0]))));

    SinCosFloat(
        x87f_store_f32(x87f_div(
            x87f_mul(x87f_neg(x87f_load_f32(angles[0])), x87f_load_f64(M_PI)),
            x87f_load_f64((double)CODUO_ANGLE_180))),
        &sinValue, &cosValue);
    axis1[0] = x87f_store_f32(
        x87f_add(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(yawTemp)),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(axis1[2]))));
    axis1[2] = x87f_store_f32(
        x87f_sub(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(axis1[2])),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(yawTemp))));
#else
    SinCosFloat((float)(((double)-angles[1] * M_PI) /
                           (double)CODUO_ANGLE_180),
                 &sinValue, &cosValue);
    yawTemp = (cosValue * axis1[0]) - (sinValue * axis1[1]);
    axis1[1] = (cosValue * axis1[1]) + (sinValue * axis1[0]);

    SinCosFloat((float)(((double)-angles[0] * M_PI) /
                           (double)CODUO_ANGLE_180),
                 &sinValue, &cosValue);
    axis1[0] = (cosValue * yawTemp) + (sinValue * axis1[2]);
    axis1[2] = (cosValue * axis1[2]) - (sinValue * yawTemp);
#endif

    float roll = (float)vectosignedpitch(axis1);
    if (axis1[1] < 0.0f) {
        if (roll < 0.0f) {
            roll += CODUO_ANGLE_180;
        } else {
            roll -= CODUO_ANGLE_180;
        }
        angles[2] = roll;
    } else {
        uint32_t rollBits;
        memcpy(&rollBits, &roll, sizeof(rollBits));
        rollBits ^= CODUO_FLOAT_SIGN_BIT;
        memcpy(&angles[2], &rollBits, sizeof(angles[2]));
    }
}

void AxisToSignedAngles(const vec3_t axis[3], vec3_t angles)
{
    vec3_t axis1;
    float sinValue;
    float cosValue;
    float yawTemp;

    vectosignedangles(axis[0], angles);

    memcpy(&axis1[0], &axis[1][0], sizeof(axis1[0]));
    memcpy(&axis1[1], &axis[1][1], sizeof(axis1[1]));
    memcpy(&axis1[2], &axis[1][2], sizeof(axis1[2]));

    /* stock multiplies by double pi then divides by double 180.0 (e.g.
     * 0x806a58d..0x806a59b); a folded pi/180 constant rounds differently */
#if EMULATE_X87
    SinCosFloat(
        x87f_store_f32(x87f_div(
            x87f_mul(x87f_neg(x87f_load_f32(angles[1])), x87f_load_f64(M_PI)),
            x87f_load_f64((double)CODUO_ANGLE_180))),
        &sinValue, &cosValue);
    yawTemp = x87f_store_f32(
        x87f_sub(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(axis1[0])),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(axis1[1]))));
    axis1[1] = x87f_store_f32(
        x87f_add(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(axis1[1])),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(axis1[0]))));

    SinCosFloat(
        x87f_store_f32(x87f_div(
            x87f_mul(x87f_neg(x87f_load_f32(angles[0])), x87f_load_f64(M_PI)),
            x87f_load_f64((double)CODUO_ANGLE_180))),
        &sinValue, &cosValue);
    axis1[0] = x87f_store_f32(
        x87f_add(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(yawTemp)),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(axis1[2]))));
    axis1[2] = x87f_store_f32(
        x87f_sub(x87f_mul(x87f_load_f32(cosValue), x87f_load_f32(axis1[2])),
                 x87f_mul(x87f_load_f32(sinValue), x87f_load_f32(yawTemp))));
#else
    SinCosFloat((float)(((double)-angles[1] * M_PI) /
                           (double)CODUO_ANGLE_180),
                 &sinValue, &cosValue);
    yawTemp = (cosValue * axis1[0]) - (sinValue * axis1[1]);
    axis1[1] = (cosValue * axis1[1]) + (sinValue * axis1[0]);

    SinCosFloat((float)(((double)-angles[0] * M_PI) /
                           (double)CODUO_ANGLE_180),
                 &sinValue, &cosValue);
    axis1[0] = (cosValue * yawTemp) + (sinValue * axis1[2]);
    axis1[2] = (cosValue * axis1[2]) - (sinValue * yawTemp);
#endif

    float roll = (float)vectosignedpitch(axis1);
    if (axis1[1] < 0.0f) {
        if (roll < 0.0f) {
            roll += CODUO_ANGLE_180;
        } else {
            roll -= CODUO_ANGLE_180;
        }
        angles[2] = roll;
    } else {
        uint32_t rollBits;
        memcpy(&rollBits, &roll, sizeof(rollBits));
        rollBits ^= CODUO_FLOAT_SIGN_BIT;
        memcpy(&angles[2], &rollBits, sizeof(angles[2]));
    }
}

qboolean PlaneFromPoints(vec4_t plane, const vec3_t point0, const vec3_t point1,
                      const vec3_t point2)
{
    vec3_t edge0;
    vec3_t edge1;

#if EMULATE_X87
    /* x87 (stock coduo_lnxded 0x0806a964), same structure as
     * CM_PlaneFromPoints: each edge
     * component is an 80-bit subtract rounded to the float slot; CrossProduct
     * and VectorNormalize are the emulated variants; plane[3] = point0 . plane
     * is three 80-bit products and two 80-bit adds rounded to float once. */
    for (int32_t axis = 0; axis < 3; ++axis) {
        edge0[axis] = x87f_store_f32(
            x87f_sub(x87f_load_f32(point1[axis]), x87f_load_f32(point0[axis])));
        edge1[axis] = x87f_store_f32(
            x87f_sub(x87f_load_f32(point2[axis]), x87f_load_f32(point0[axis])));
    }

    CrossProduct(edge1, edge0, plane);
    if (VectorNormalize(plane) == 0.0) {
        return qfalse;
    }

    plane[3] = x87f_store_f32(x87f_add(
        x87f_add(x87f_mul(x87f_load_f32(point0[0]), x87f_load_f32(plane[0])),
                 x87f_mul(x87f_load_f32(point0[1]), x87f_load_f32(plane[1]))),
        x87f_mul(x87f_load_f32(point0[2]), x87f_load_f32(plane[2]))));
    return qtrue;
#else
    edge0[0] = point1[0] - point0[0];
    edge0[1] = point1[1] - point0[1];
    edge0[2] = point1[2] - point0[2];
    edge1[0] = point2[0] - point0[0];
    edge1[1] = point2[1] - point0[1];
    edge1[2] = point2[2] - point0[2];

    CrossProduct(edge1, edge0, plane);
    if (VectorNormalize(plane) == 0.0) {
        return qfalse;
    }

    plane[3] = (point0[0] * plane[0]) + (point0[1] * plane[1]) +
               (point0[2] * plane[2]);
    return qtrue;
#endif
}

void ProjectPointOnPlane(vec3_t dst, const vec3_t point, const vec3_t normal)
{
    /* ORIGINAL_BINARY_BUG: stock scales both the normal dot product and the
     * normal itself by 1/|normal|^2, producing a |normal|^-4 projection for
     * non-unit normals. Preserved; see
     * original-bugs/game-project-point-on-plane-double-inverse.md. */
#if EMULATE_X87
    /* stock 0x806aa4e: every named local is one 80-bit x87 chain rounded once
     * at its float slot. lengthSq = n·n (fstps -0x30); invLengthSq = 1/lengthSq
     * (fld1;fdivs;fstps -0x2c); projected = (n·point)*invLengthSq — the dot is
     * kept 80-bit and multiplied by invLengthSq before the single store (fstps
     * -0xc); each scaled[i] = n[i]*invLengthSq is rounded to its own slot
     * (0x806aac1..0x806aae5); dst[i] = point[i] - projected*scaled[i]
     * (fmuls;fsubp;fstps). */
    const float lengthSq = x87f_store_f32(
        x87f_add(x87f_add(x87f_mul(x87f_load_f32(normal[0]),
                                   x87f_load_f32(normal[0])),
                          x87f_mul(x87f_load_f32(normal[1]),
                                   x87f_load_f32(normal[1]))),
                 x87f_mul(x87f_load_f32(normal[2]), x87f_load_f32(normal[2]))));
    const float invLengthSq = x87f_store_f32(
        x87f_div(x87f_load_f32(1.0f), x87f_load_f32(lengthSq)));
    const float projected = x87f_store_f32(x87f_mul(
        x87f_add(x87f_add(x87f_mul(x87f_load_f32(normal[0]),
                                   x87f_load_f32(point[0])),
                          x87f_mul(x87f_load_f32(normal[1]),
                                   x87f_load_f32(point[1]))),
                 x87f_mul(x87f_load_f32(normal[2]), x87f_load_f32(point[2]))),
        x87f_load_f32(invLengthSq)));
    const float scaled0 = x87f_store_f32(
        x87f_mul(x87f_load_f32(normal[0]), x87f_load_f32(invLengthSq)));
    const float scaled1 = x87f_store_f32(
        x87f_mul(x87f_load_f32(normal[1]), x87f_load_f32(invLengthSq)));
    const float scaled2 = x87f_store_f32(
        x87f_mul(x87f_load_f32(normal[2]), x87f_load_f32(invLengthSq)));

    dst[0] = x87f_store_f32(x87f_sub(x87f_load_f32(point[0]),
                                     x87f_mul(x87f_load_f32(projected),
                                              x87f_load_f32(scaled0))));
    dst[1] = x87f_store_f32(x87f_sub(x87f_load_f32(point[1]),
                                     x87f_mul(x87f_load_f32(projected),
                                              x87f_load_f32(scaled1))));
    dst[2] = x87f_store_f32(x87f_sub(x87f_load_f32(point[2]),
                                     x87f_mul(x87f_load_f32(projected),
                                              x87f_load_f32(scaled2))));
#else
    /* stock 0x806aa82 rounds the squared length to a float slot before the
     * reciprocal */
    float lengthSq = (normal[0] * normal[0]) + (normal[1] * normal[1]) +
                     (normal[2] * normal[2]);
    float invLengthSq = 1.0f / lengthSq;
    float projected =
        ((normal[0] * point[0]) + (normal[1] * point[1]) +
         (normal[2] * point[2])) *
        invLengthSq;
    /* stock 0x806aac1..0x806aae5 rounds each normal[i]*invLengthSq factor to
     * a float slot before the dst subtractions */
    float scaled0 = normal[0] * invLengthSq;
    float scaled1 = normal[1] * invLengthSq;
    float scaled2 = normal[2] * invLengthSq;

    dst[0] = point[0] - (projected * scaled0);
    dst[1] = point[1] - (projected * scaled1);
    dst[2] = point[2] - (projected * scaled2);
#endif
}

qboolean BoxDistSqrdExceeds(const vec3_t point0, const vec3_t point1,
                      const vec3_t point2, float epsilon)
{
    vec3_t delta0;
    vec3_t delta1;
    float sum = 0.0f;

    for (int32_t axis = 0; axis < 3; ++axis) {
        delta0[axis] = point0[axis] - point2[axis];
        delta1[axis] = point1[axis] - point2[axis];
        if (!(0.0f >= (delta0[axis] * delta1[axis]))) {
            float delta0Sq = delta0[axis] * delta0[axis];
            float delta1Sq = delta1[axis] * delta1[axis];
            sum += delta0Sq > delta1Sq ? delta1Sq : delta0Sq;
        }
    }

    return epsilon < sum ? qtrue : qfalse;
}

void NormalToLatLong(const vec3_t normal, uint8_t encoded[2])
{
    if (normal[0] == 0.0f && normal[1] == 0.0f) {
        encoded[0] = normal[2] > 0.0f ? 0 : CODUO_NORMAL_BYTE_NEGATIVE_Z;
        encoded[1] = 0;
        return;
    }

    /* stock 0x806ad01..0x806ad17 (and the acos twin): the radians value is
     * multiplied by double 180.0, divided by double pi, then multiplied by
     * the float-folded 255.0f/360.0f (rodata 0x80dd8f0), NOT by 256/360 */
#if EMULATE_X87
    /* x87-faithful transcription (stock coduo_lnxded 0x0806ac7e): the libm
     * radians (atan2 / acos) are scaled by 180.0 then divided by pi then
     * multiplied by the float-folded 255/360, all in 80-bit, and truncated toward
     * zero to int (the C int16 cast). atan2/acos are libm calls left native. */
    uint8_t lng = (uint8_t)(int16_t)x87f_store_i32_trunc(x87f_mul(
        x87f_div(x87f_mul(x87f_load_f64(atan2(normal[1], normal[0])),
                          x87f_load_f64((double)CODUO_ANGLE_180)),
                 x87f_load_f64(M_PI)),
        x87f_load_f32((float)CODUO_NORMAL_ENCODE_SCALE /
                      (float)CODUO_NORMAL_BYTE_SCALE)));
    uint8_t lat = (uint8_t)(int16_t)x87f_store_i32_trunc(x87f_mul(
        x87f_div(x87f_mul(x87f_load_f64(acos(normal[2])),
                          x87f_load_f64((double)CODUO_ANGLE_180)),
                 x87f_load_f64(M_PI)),
        x87f_load_f32((float)CODUO_NORMAL_ENCODE_SCALE /
                      (float)CODUO_NORMAL_BYTE_SCALE)));
#elif defined(__x86_64__)
    static const double degrees = (double)CODUO_ANGLE_180;
    static const double pi = M_PI;
    static const double scale =
        (double)((float)CODUO_NORMAL_ENCODE_SCALE /
                 (float)CODUO_NORMAL_BYTE_SCALE);
    double lngRadians = atan2(normal[1], normal[0]);
    uint8_t lng = (uint8_t)CODUO_X87_TRUNCATE_I16(
        ((long double)lngRadians * (long double)degrees /
          (long double)pi) *
        (long double)scale);
    double latRadians = acos(normal[2]);
    uint8_t lat = (uint8_t)CODUO_X87_TRUNCATE_I16(
        ((long double)latRadians * (long double)degrees /
          (long double)pi) *
        (long double)scale);
#else
    uint8_t lng = (uint8_t)(int16_t)(((atan2(normal[1], normal[0]) *
                                       (double)CODUO_ANGLE_180) /
                                      M_PI) *
                                     ((float)CODUO_NORMAL_ENCODE_SCALE /
                                      (float)CODUO_NORMAL_BYTE_SCALE));
    uint8_t lat = (uint8_t)(int16_t)(((acos(normal[2]) *
                                       (double)CODUO_ANGLE_180) /
                                      M_PI) *
                                     ((float)CODUO_NORMAL_ENCODE_SCALE /
                                      (float)CODUO_NORMAL_BYTE_SCALE));
#endif

    encoded[0] = lat;
    encoded[1] = lng;
}

void Vec10Copy(const float in[10], float out[10])
{
    for (int32_t index = 0; index < 10; ++index) {
        memcpy(&out[index], &in[index], sizeof(out[index]));
    }
}

long double RoundFloat(float value)
{
#if EMULATE_X87
    /* value + 0.5 in 80-bit rounded to double64 for the CRT floor, floor result
     * narrowed to float (stock coduo_lnxded 0x0806ae2f); floor is exact. */
    return (float)floor(
        x87f_store_f64(x87f_add(x87f_load_f32(value), x87f_load_f32(0.5f))));
#else
    return (float)floor((double)value + 0.5);
#endif
}

long double ColorNormalize(const vec3_t source, vec3_t dest)
{
    float max;
    memcpy(&max, &source[0], sizeof(max));

    if (max < source[1]) {
        memcpy(&max, &source[1], sizeof(max));
    }
    if (max < source[2]) {
        memcpy(&max, &source[2], sizeof(max));
    }

    if (max == 0.0f) {
        const uint32_t oneBits = UINT32_C(0x3f800000); /* 1.0f */
        memcpy(&dest[2], &oneBits, sizeof(dest[2]));
        memcpy(&dest[1], &oneBits, sizeof(dest[1]));
        memcpy(&dest[0], &oneBits, sizeof(dest[0]));
        return 0.0f;
    }

#if EMULATE_X87
    /* reciprocal then three 80-bit scales, each rounded to float; the
     * max-selection compares are plain float compares. */
    float scale = x87f_store_f32(x87f_div(x87f_load_f32(1.0f), x87f_load_f32(max)));
    dest[0] = x87f_store_f32(x87f_mul(x87f_load_f32(source[0]),
                                      x87f_load_f32(scale)));
    dest[1] = x87f_store_f32(x87f_mul(x87f_load_f32(source[1]),
                                      x87f_load_f32(scale)));
    dest[2] = x87f_store_f32(x87f_mul(x87f_load_f32(source[2]),
                                      x87f_load_f32(scale)));
#else
    float scale = 1.0f / max;
    dest[0] = source[0] * scale;
    dest[1] = source[1] * scale;
    dest[2] = source[2] * scale;
#endif
    return max;
}

void VectorRotateAngles(const vec3_t point, const vec3_t angles, vec3_t out)
{
    vec3_t previous;
    vec3_t rotated;
    int32_t axisPairs[6] = {1, 2, 2, 0, 0, 1};

    memcpy(&previous[0], &point[0], sizeof(previous[0]));
    memcpy(&previous[1], &point[1], sizeof(previous[1]));
    memcpy(&previous[2], &point[2], sizeof(previous[2]));
    memcpy(&rotated[0], &previous[0], sizeof(rotated[0]));
    memcpy(&rotated[1], &previous[1], sizeof(rotated[1]));
    memcpy(&rotated[2], &previous[2], sizeof(rotated[2]));

    for (int32_t angleIndex = 0; angleIndex < 3; ++angleIndex) {
        if (angles[angleIndex] != 0.0f) {
            double sinValue;
            double cosValue;
#if EMULATE_X87
            /* (angles[i]*pi)/180 in 80-bit rounded to double for SinCosDouble;
             * each rotated component multiplies the float lanes by the full
             * double sin/cos in 80-bit and rounds to float. */
            double angle = x87f_store_f64(
                x87f_div(x87f_mul(x87f_load_f32(angles[angleIndex]),
                                  x87f_load_f64(M_PI)),
                         x87f_load_f64((double)CODUO_ANGLE_180)));
#else
            double angle = ((double)angles[angleIndex] * M_PI) /
                           (double)CODUO_ANGLE_180;
#endif

            SinCosDouble(angle, &sinValue, &cosValue);

            int32_t first = axisPairs[angleIndex * 2];
            int32_t second = axisPairs[(angleIndex * 2) + 1];
            /* stock 0x806b010/0x806b01e multiplies the float components by
             * the full double sin/cos (fmul QWORD) with no float rounding of
             * the sin/cos values */
#if EMULATE_X87
            rotated[first] = x87f_store_f32(
                x87f_sub(x87f_mul(x87f_load_f32(previous[first]),
                                  x87f_load_f64(cosValue)),
                         x87f_mul(x87f_load_f32(previous[second]),
                                  x87f_load_f64(sinValue))));
            rotated[second] = x87f_store_f32(
                x87f_add(x87f_mul(x87f_load_f32(previous[second]),
                                  x87f_load_f64(cosValue)),
                         x87f_mul(x87f_load_f32(previous[first]),
                                  x87f_load_f64(sinValue))));
#else
            rotated[first] = (previous[first] * cosValue) -
                             (previous[second] * sinValue);
            rotated[second] = (previous[second] * cosValue) +
                              (previous[first] * sinValue);
#endif
        }

        memcpy(&previous[0], &rotated[0], sizeof(previous[0]));
        memcpy(&previous[1], &rotated[1], sizeof(previous[1]));
        memcpy(&previous[2], &rotated[2], sizeof(previous[2]));
    }

    memcpy(&out[0], &rotated[0], sizeof(out[0]));
    memcpy(&out[1], &rotated[1], sizeof(out[1]));
    memcpy(&out[2], &rotated[2], sizeof(out[2]));
}

void VectorRotateAnglesAroundPoint(const vec3_t point, const vec3_t angles, const vec3_t origin,
                  vec3_t out)
{
    vec3_t offset;
    vec3_t rotated;

    offset[0] = point[0] - origin[0];
    offset[1] = point[1] - origin[1];
    offset[2] = point[2] - origin[2];
    VectorRotateAngles(offset, angles, rotated);

    out[0] = rotated[0] + origin[0];
    out[1] = rotated[1] + origin[1];
    out[2] = rotated[2] + origin[2];
}

void VectorPolar(vec3_t out, float radius, float angle)
{
    float sin0;
    float cos0;
    float sin1;
    float cos1;

    SinCosFloat(angle, &sin0, &cos0);
    SinCosFloat(angle, &sin1, &cos1);

#if EMULATE_X87
    /* each component is an 80-bit product chain rounded to float (stock
     * coduo_lnxded 0x0806b113). */
    out[0] = x87f_store_f32(x87f_mul(
        x87f_mul(x87f_load_f32(radius), x87f_load_f32(cos0)),
        x87f_load_f32(cos1)));
    out[1] = x87f_store_f32(x87f_mul(
        x87f_mul(x87f_load_f32(radius), x87f_load_f32(sin0)),
        x87f_load_f32(cos1)));
    out[2] =
        x87f_store_f32(x87f_mul(x87f_load_f32(radius), x87f_load_f32(sin1)));
#else
    out[0] = radius * cos0 * cos1;
    out[1] = radius * sin0 * cos1;
    out[2] = radius * sin1;
#endif
}

void VectorSnap(vec3_t vector)
{
#if EMULATE_X87
    /* each component + 0.5 in 80-bit rounded to double64 for floor, floor result
     * narrowed to float (stock coduo_lnxded 0x0806b17a). */
    vector[0] = (float)floor(
        x87f_store_f64(x87f_add(x87f_load_f32(vector[0]), x87f_load_f32(0.5f))));
    vector[1] = (float)floor(
        x87f_store_f64(x87f_add(x87f_load_f32(vector[1]), x87f_load_f32(0.5f))));
    vector[2] = (float)floor(
        x87f_store_f64(x87f_add(x87f_load_f32(vector[2]), x87f_load_f32(0.5f))));
#else
    vector[0] = (float)floor((double)vector[0] + 0.5);
    vector[1] = (float)floor((double)vector[1] + 0.5);
    vector[2] = (float)floor((double)vector[2] + 0.5);
#endif
}

void _Vector5Add(const float left[5], const float right[5], float out[5])
{
    for (int32_t index = 0; index < 5; ++index) {
        out[index] = left[index] + right[index];
    }
}

void _Vector5Scale(const float source[5], float scale, float out[5])
{
    for (int32_t index = 0; index < 5; ++index) {
        out[index] = source[index] * scale;
    }
}

void _Vector53Copy(const vec3_t source, vec3_t dest)
{
    memcpy(&dest[0], &source[0], sizeof(dest[0]));
    memcpy(&dest[1], &source[1], sizeof(dest[1]));
    memcpy(&dest[2], &source[2], sizeof(dest[2]));
}

long double round_decimal(float value, int32_t decimals)
{
#if EMULATE_X87
    /* x87-faithful transcription (stock coduo_lnxded 0x0806b2e2): each pow return
     * is multiplied in 80-bit (fmul m64) against (double)value / whole with no
     * intermediate double rounding of the scale; scaledValue rounds to double64,
     * the final product rounds to float. modf and the whole+-1 adjustments are
     * exact double ops. pow is a libm call left native (arm64-vs-x86 libm
     * concern, same as the other transcendentals). */
    double scaledValue = x87f_store_f64(x87f_mul(
        x87f_load_f64(pow(10.0, (double)decimals)), x87f_load_f64((double)value)));
    double whole;
    double frac = modf(scaledValue, &whole);

    if (frac >= 0.5) {
        whole += 1.0;
    } else if (frac <= -0.5) {
        whole -= 1.0;
    }

    return (long double)x87f_store_f32(x87f_mul(
        x87f_load_f64(pow(0.1, (double)decimals)), x87f_load_f64(whole)));
#else
    /* stock 0x806b309 multiplies the raw x87 pow return (no intermediate
     * double rounding of the scale) */
    double scaledValue = pow(10.0, (double)decimals) * (double)value;
    double whole;
    double frac = modf(scaledValue, &whole);

    if (frac >= 0.5) {
        whole += 1.0;
    } else if (frac <= -0.5) {
        whole -= 1.0;
    }

    return (float)(pow(0.1, (double)decimals) * whole);
#endif
}

long double PitchForYawOnNormal(float angleDegrees, const vec3_t normal)
{
    vec3_t forward;
    vec3_t projected;

    YawVectors(angleDegrees, forward, NULL);
    ProjectPointOnPlane(projected, forward, normal);
    return vectopitch(projected);
}

void Rand_Init(uint32_t seed)
{
    randSeed = seed;
}

long double flrand(float min, float max)
{
    randSeed = (randSeed * 214013U) + 2531011U;
#if EMULATE_X87
    /* x87-faithful transcription (stock coduo_lnxded 0x0806b3ea): (randSeed>>17)
     * is fild'd and rounded to a float slot; then ((max-min)*that)/32768 + min is
     * evaluated left-to-right in 80-bit and rounded to float. */
    float randFloat = x87f_store_f32(x87f_load_i32((int32_t)(randSeed >> 17)));
    float result = x87f_store_f32(x87f_add(
        x87f_div(x87f_mul(x87f_sub(x87f_load_f32(max), x87f_load_f32(min)),
                          x87f_load_f32(randFloat)),
                 x87f_load_f32(32768.0f)),
        x87f_load_f32(min)));
    return result;
#else
    /* stock 0x806b432 rounds the result to float before the x87 return */
    float result = (((max - min) * (float)(randSeed >> 17)) / 32768.0f) + min;
    return result;
#endif
}

int32_t irand(int32_t min, int32_t max)
{
    randSeed = (randSeed * 214013U) + 2531011U;
    int32_t range = max - min;
    int32_t random = (int32_t)(randSeed >> 17);
    int32_t product = range * random;
    int32_t scaled = coduo_int32_from_bits(
        coduo_int32_sar_bits((uint32_t)product, 15U));

    return scaled + min;
}

long double Q_SwayRand(float frequency0, float frequency1, float time)
{
#if EMULATE_X87
    /* x87-faithful transcription (stock coduo_lnxded 0x0806b480): each angle
     * frequency*(time/1000)*pi is an 80-bit chain rounded to double for the libm
     * sin/cos; the (angle+angle) doublings are single-op double adds; the final
     * cos*sinValue is 80-bit rounded to float. sin/cos are libm calls left
     * native. */
    double angle = x87f_store_f64(x87f_mul(
        x87f_mul(x87f_load_f32(frequency0),
                 x87f_div(x87f_load_f32(time), x87f_load_f64(1000.0))),
        x87f_load_f64(M_PI)));
    double sinValue = sin(angle + angle);

    angle = x87f_store_f64(x87f_mul(
        x87f_mul(x87f_load_f32(frequency1),
                 x87f_div(x87f_load_f32(time), x87f_load_f64(1000.0))),
        x87f_load_f64(M_PI)));
    return (long double)x87f_store_f32(x87f_mul(
        x87f_load_f64(cos(angle + angle)), x87f_load_f64(sinValue)));
#else
    double angle = ((double)frequency0 * ((double)time / 1000.0)) * M_PI;
    double sinValue = sin(angle + angle);

    angle = ((double)frequency1 * ((double)time / 1000.0)) * M_PI;
    /* stock 0x806b4cd multiplies the raw x87 cos return by sinValue (no
     * intermediate double rounding of the cosine) */
    return (float)((long double)cos(angle + angle) * (long double)sinValue);
#endif
}

void SinCosFloat(float angle, float *sinValue, float *cosValue)
{
#if defined(__i386__) || defined(__x86_64__)
    coduo_x87_sincosf(angle, sinValue, cosValue);
#elif defined(CODUO_ENGINE_ALLOW_NON_X87_FLOAT)
    *cosValue = cosf(angle);
    *sinValue = sinf(angle);
#else
#error "SinCosFloat requires x87 fsincos or CODUO_ENGINE_ALLOW_NON_X87_FLOAT"
#endif
}

void SinCosDouble(double angle, double *sinValue, double *cosValue)
{
    *sinValue = sin(angle);
    *cosValue = cos(angle);
}

const char *Com_StringContains(const char *haystack, const char *needle,
                         qboolean caseSensitive)
{
    int32_t maxOffset =
        (int32_t)strlen(haystack) - (int32_t)strlen(needle);

    for (int32_t offset = 0; offset <= maxOffset;
         ++offset, ++haystack) {
        int32_t index;
        for (index = 0; needle[index] != '\0'; ++index) {
            if (caseSensitive == qfalse) {
                if (toupper(coduo_ctype_signed_byte_arg(haystack[index])) !=
                    toupper(coduo_ctype_signed_byte_arg(needle[index]))) {
                    break;
                }
            } else if (haystack[index] != needle[index]) {
                break;
            }
        }

        if (needle[index] == '\0') {
            return haystack;
        }
    }

    return NULL;
}

qboolean Com_Filter(const char *filter, const char *name,
                      qboolean caseSensitive)
{
    /* ORIGINAL_BINARY_BUG: the post-'*' literal run has no 1,024-byte bound,
     * '?' advances past a short candidate without a NUL test, and a matching
     * unterminated class advances beyond the filter terminator. All three
     * retail behaviors are preserved; see the coduomp-com-filter-* reports in
     * original-bugs/. */
    char chunk[CODUO_WILDCARD_CHUNK_SIZE];

    for (;;) {
        if (*filter == '\0') {
            return qtrue;
        }

        if (*filter == '*') {
            int32_t chunkLength = 0;
            while (*++filter != '\0' && *filter != '*' && *filter != '?') {
                chunk[chunkLength] = *filter;
                ++chunkLength;
            }
            chunk[chunkLength] = '\0';

            if (chunk[0] != '\0') {
                const char *match =
                    Com_StringContains(name, chunk, caseSensitive);
                if (match == NULL) {
                    return qfalse;
                }
                name = match + strlen(chunk);
            }
            continue;
        }

        if (*filter == '?') {
            ++filter;
            ++name;
            continue;
        }

        if (*filter == '[' && filter[1] == '[') {
            ++filter;
            continue;
        }

        if (*filter == '[') {
            qboolean matched = qfalse;
            ++filter;

            while (*filter != '\0' && matched == qfalse &&
                   (*filter != ']' || filter[1] == ']')) {
                if (filter[1] == '-' && filter[2] != '\0' &&
                    (filter[2] != ']' || filter[3] == ']')) {
                    if (caseSensitive != qfalse) {
                        if (*filter <= *name && *name <= filter[2]) {
                            matched = qtrue;
                        }
                    } else {
                        int32_t nameChar =
                            toupper(coduo_ctype_signed_byte_arg(*name));
                        if (toupper(coduo_ctype_signed_byte_arg(*filter)) <=
                                nameChar &&
                            nameChar <= toupper(coduo_ctype_signed_byte_arg(
                                            filter[2]))) {
                            matched = qtrue;
                        }
                    }
                    filter += 3;
                } else {
                    if (caseSensitive != qfalse) {
                        if (*filter == *name) {
                            matched = qtrue;
                        }
                    } else if (toupper(coduo_ctype_signed_byte_arg(*filter)) ==
                               toupper(coduo_ctype_signed_byte_arg(*name))) {
                        matched = qtrue;
                    }
                    ++filter;
                }
            }

            if (matched == qfalse) {
                return qfalse;
            }

            while (*filter != '\0' && (*filter != ']' || filter[1] == ']')) {
                ++filter;
            }
            ++filter;
            ++name;
            continue;
        }

        if (caseSensitive == qfalse) {
            if (toupper(coduo_ctype_signed_byte_arg(*filter)) !=
                toupper(coduo_ctype_signed_byte_arg(*name))) {
                return qfalse;
            }
        } else if (*filter != *name) {
            return qfalse;
        }

        ++filter;
        ++name;
    }
}

qboolean Com_FilterPath(const char *filter, const char *name,
                        qboolean caseSensitive)
{
    char normalizedFilter[CODUO_PATH_MATCH_BUFFER_SIZE];
    char normalizedName[CODUO_PATH_MATCH_BUFFER_SIZE];
    int32_t index;

    for (index = 0; index < CODUO_PATH_MATCH_BUFFER_SIZE - 1 &&
                    filter[index] != '\0';
         ++index) {
        normalizedFilter[index] =
            (filter[index] == '\\' || filter[index] == ':') ? '/'
                                                            : filter[index];
    }
    normalizedFilter[index] = '\0';

    for (index = 0;
         index < CODUO_PATH_MATCH_BUFFER_SIZE - 1 && name[index] != '\0';
         ++index) {
        normalizedName[index] =
            (name[index] == '\\' || name[index] == ':') ? '/' : name[index];
    }
    normalizedName[index] = '\0';

    return Com_Filter(normalizedFilter, normalizedName, caseSensitive);
}

uint32_t Com_HashKey(const char *text, int32_t length)
{
    int32_t hash = 0;

    for (int32_t index = 0; index < length && text[index] != '\0'; ++index) {
        hash += (index + CODUO_HASH_POSITION_BIAS) * text[index];
    }

    return (uint32_t)(hash ^ (hash >> 10) ^ (hash >> 20));
}

time_t Com_RealTime(qtime_t *qtime)
{
    time_t now = time(NULL);

    if (qtime != NULL) {
        struct tm *localTime = localtime(&now);
        if (localTime != NULL) {
            qtime->tm_sec = localTime->tm_sec;
            qtime->tm_min = localTime->tm_min;
            qtime->tm_hour = localTime->tm_hour;
            qtime->tm_mday = localTime->tm_mday;
            qtime->tm_mon = localTime->tm_mon;
            qtime->tm_year = localTime->tm_year;
            qtime->tm_wday = localTime->tm_wday;
            qtime->tm_yday = localTime->tm_yday;
            qtime->tm_isdst = localTime->tm_isdst;
        }
    }

    return now;
}

int32_t Q_rint(float value)
{
#if EMULATE_X87
    /* value + 0.5f in 80-bit truncated toward zero (stock fadds; fistpl under a
     * 0x0c00 truncate control word), reproduced with the shim's truncating
     * store — the arm64 path where the native x87 inline asm is unavailable. */
    return x87f_store_i32_trunc(
        x87f_add(x87f_load_f32(value), x87f_load_f32(0.5f)));
#elif CODUO_ENGINE_HAS_X87_INLINE_ASM
    const float half = 0.5f;
    uint16_t savedControlWord;
    uint16_t truncateControlWord;
    int32_t result;

    __asm__ __volatile__("fnstcw %0"
                         : "=m"(savedControlWord)
                         :
                         : "memory");
    truncateControlWord = (uint16_t)(savedControlWord | 0x0c00);

    __asm__ __volatile__(
        "flds %1\n\t"
        "fadds %2\n\t"
        "fldcw %3\n\t"
        "fistpl %0\n\t"
        "fldcw %4"
        : "=m"(result)
        : "m"(value), "m"(half), "m"(truncateControlWord),
          "m"(savedControlWord)
        : "memory");

    return result;
#elif defined(CODUO_ENGINE_ALLOW_NON_X87_FLOAT)
    return (int32_t)(value + 0.5f);
#else
#error "Q_rint requires x87 inline assembly for original integer conversion"
#endif
}
