#include <math.h>
#include <stddef.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <strings.h>

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

#include "../core_math/core_math_private.h"
#include "../core_memory/core_memory_private.h"
#include "../core_runtime/core_runtime_private.h"
#include "../model_assets/xmodel_asset_private.h"
#include "coduo_engine_structs.h"
#include "cm_collision_globals_private.h"
#include "cm_terrain_private.h"
#include "cm_world_sector_private.h"

#define CODUO_CM_WINDING_CLIP_EPSILON 0.1f
#define CODUO_CM_WINDING_MIN_AREA 1.0f
#define CODUO_CM_WINDING_AXIS_NORMAL 1.0f

enum {
    CM_PATCH_VECTOR_AXIS_COUNT = 3,
    CODUO_CM_PATCH_POINT_GRID_SIZE = 129,
    CODUO_CM_PATCH_PLANE_GRID_SIZE = 129,
    CODUO_CM_PATCH_EDGE_PLANE_PUSH_DISTANCE = 4,
    CODUO_CM_PATCH_PLANE_SIDE_FRONT = 0,
    CODUO_CM_PATCH_PLANE_SIDE_BACK = 1,
    CODUO_CM_PATCH_PLANE_SIDE_ON = 2,
    CODUO_CM_PATCH_PLANE_SIDE_CROSS = -2,
    CODUO_CM_WINDING_MAX_POINTS = 64,
    CODUO_CM_WINDING_TEMP_POINT_LIMIT = 65,
    CODUO_CM_WINDING_MERGE_POINT_LIMIT = 128,
    CODUO_CM_WINDING_SPLIT_EXTRA_POINTS = 4,
    CODUO_CM_WINDING_FREED_SENTINEL = -0x21522153,
    CODUO_CM_WINDING_WORLD_COORD_LIMIT = 131072,
    CODUO_CM_PATCH_EDGE_BOTTOM = 0,
    CODUO_CM_PATCH_EDGE_RIGHT = 1,
    CODUO_CM_PATCH_EDGE_TOP = 2,
    CODUO_CM_PATCH_EDGE_LEFT = 3,
    CODUO_CM_PATCH_EDGE_DIAGONAL_DESCENDING = 4,
    CODUO_CM_PATCH_EDGE_DIAGONAL_ASCENDING = 5,
    CODUO_CM_PATCH_BORDER_INWARD_LOWER_TRIANGLE = 0,
    CODUO_CM_PATCH_BORDER_INWARD_UPPER_TRIANGLE = 1,
    CODUO_CM_PATCH_BORDER_INWARD_QUAD = -1,
    CODUO_CM_WINDING_MIN_POINT_COUNT = 3,
    CODUO_CM_PATCH_MAX_FACETS = 1024,
    CODUO_CM_PATCH_BEVEL_AXIS_MIN = -1,
    CODUO_CM_PATCH_BEVEL_AXIS_MAX = 1,
    CODUO_CM_MISC_MODEL_KEY_BUFFER_SIZE = 64,
    CODUO_CM_MISC_MODEL_VALUE_BUFFER_SIZE = 64,
    CODUO_CM_MISC_MODEL_NAME_BUFFER_SIZE = 64,
    CODUO_CM_ENTITYNUM_WORLD = 1022,
    CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT = 3,
    CODUO_CM_MODEL_EDGE_LIMIT = 65536,
    CODUO_CM_MODEL_VERTEX_REMAP_LIMIT = 65536,
};

#define CODUO_CM_MODEL_COPLANAR_EDGE_EPSILON 6.4e-09f
#define CODUO_CM_MODEL_TRIANGLE_PLANE_NEGATIVE_EPSILON -0.001f
#define CODUO_CM_MODEL_TRIANGLE_PLANE_POSITIVE_EPSILON 0.001f
#define CODUO_CM_PATCH_WRAP_EPSILON 0.1
#define CODUO_CM_PATCH_POINT_COMPARE_EPSILON 0.1
#define CODUO_CM_PATCH_PLANE_NORMAL_EPSILON 0.0001
#define CODUO_CM_PATCH_SNAP_VECTOR_EPSILON 0.0001
#define CODUO_CM_PATCH_PLANE_DIST_EPSILON 0.02
#define CODUO_CM_PATCH_POINT_PLANE_EPSILON 0.1
#define CODUO_CM_WINDING_COLINEAR_DOT_EPSILON 0.999
/* These are double constants in the stock binary (fld QWORD), not float. It
 * matters for BEVEL_PLANE_EPSILON: the bevel `dist > EPSILON` / `dist < -EPSILON`
 * classification compares `dist` against double 0.1 (0x3fb9999999999a), whereas
 * a float 0.1f (0x3dcccccd ~ 0.100000001) is a different value and flips the
 * branch at the boundary (e.g. dist == 0.1f: `> 0.1f` is false but `> 0.1` is
 * true). The three ChopWindingInPlace(...,float epsilon) call sites coerce it
 * back to 0.1f, so those are unaffected. 0.5 is exact in both widths, so
 * BEVEL_EDGE_LENGTH_MIN is unchanged in value; kept as double to match the
 * stock constant. */
#define CODUO_CM_PATCH_BEVEL_PLANE_EPSILON 0.1
#define CODUO_CM_PATCH_BEVEL_EDGE_LENGTH_MIN 0.5

typedef struct coduo_cm_triangle_edge_s {
    int16_t first;
    int16_t second;
} coduo_cm_triangle_edge_t;

typedef struct coduo_cm_collision_vertex_s {
    int32_t checkcount;
    vec3_t position;
} coduo_cm_collision_vertex_t;

typedef struct coduo_cm_collision_edge_s {
    int32_t checkcount;
    vec3_t origin;
    vec3_t radialAxes[2];
    vec3_t unitDirection;
    float length;
} coduo_cm_collision_edge_t;

typedef struct coduo_cm_collision_triangle_s {
    vec4_t plane;
    vec4_t svec;
    vec4_t tvec;
    coduo_cm_collision_vertex_t
        *vertices[CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT];
    coduo_cm_collision_edge_t
        *oppositeEdges[CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT];
} coduo_cm_collision_triangle_t;

typedef struct coduo_cm_collision_triangle_soup_s {
    uint16_t triangleCount;
    uint8_t negativeZOnly;
    uint8_t padding03;
    coduo_cm_collision_triangle_t triangles[];
} coduo_cm_collision_triangle_soup_t;

#if UINTPTR_MAX == UINT32_MAX
_Static_assert(sizeof(coduo_cm_triangle_edge_t) == 0x04,
               "coduo_cm_triangle_edge_t size mismatch");
_Static_assert(sizeof(coduo_cm_collision_vertex_t) == 0x10,
               "coduo_cm_collision_vertex_t size mismatch");
_Static_assert(sizeof(coduo_cm_collision_edge_t) == 0x38,
               "coduo_cm_collision_edge_t size mismatch");
_Static_assert(offsetof(coduo_cm_collision_triangle_t, vertices) == 0x30,
               "coduo_cm_collision_triangle_t.vertices offset mismatch");
_Static_assert(offsetof(coduo_cm_collision_triangle_t, oppositeEdges) ==
                   0x30 + 3 *
                              sizeof(coduo_cm_collision_vertex_t *),
               "coduo_cm_collision_triangle_t.oppositeEdges offset mismatch");
_Static_assert(sizeof(coduo_cm_collision_triangle_t) == 0x48,
               "coduo_cm_collision_triangle_t size mismatch");
_Static_assert(offsetof(coduo_cm_collision_triangle_soup_t, triangles) ==
                   ((0x04 + _Alignof(coduo_cm_collision_triangle_t) - 1) /
                    _Alignof(coduo_cm_collision_triangle_t)) *
                       _Alignof(coduo_cm_collision_triangle_t),
               "coduo_cm_collision_triangle_soup_t.triangles offset mismatch");
#endif

typedef struct cGrid_s {
    int32_t width;
    int32_t height;
    qboolean wrapWidth;
    qboolean wrapHeight;
    vec3_t points[CODUO_CM_PATCH_POINT_GRID_SIZE]
                 [CODUO_CM_PATCH_POINT_GRID_SIZE];
} cGrid_t;

typedef struct winding_s {
    int32_t numpoints;
    vec3_t p[];
} winding_t;

#if UINTPTR_MAX == UINT32_MAX
_Static_assert(offsetof(cGrid_t, points) == 0x10,
               "cGrid_t.points offset mismatch");
_Static_assert(sizeof(((cGrid_t *)0)->points[0]) ==
                   0x60c,
               "cGrid_t.points column stride mismatch");
_Static_assert(offsetof(winding_t, numpoints) == 0x00,
               "winding_t.numpoints offset mismatch");
_Static_assert(offsetof(winding_t, p) == 0x04,
               "winding_t.p offset mismatch");
_Static_assert(sizeof(winding_t) == 0x04,
               "i386 winding_t header size mismatch");
#endif

static patchPlane_t cm_patchPlanes[CM_PATCH_MAX_PLANES];
static int32_t cm_patchPlaneCount;
static int32_t cm_windingActiveCount;
static int32_t cm_windingPeakActiveCount;

static float cm_windingSplitDist;
static float cm_windingChopDist;
static int32_t cm_patchPlaneGrid[CODUO_CM_PATCH_PLANE_GRID_SIZE]
                                [CODUO_CM_PATCH_PLANE_GRID_SIZE][2];
static cGrid_t cm_patchWorkGrid;
static facet_t cm_patchFacets[CODUO_CM_PATCH_MAX_FACETS];
static coduo_cm_triangle_edge_t cm_modelEdges[CODUO_CM_MODEL_EDGE_LIMIT];
static int16_t cm_modelEdgeOppositeVertices[CODUO_CM_MODEL_EDGE_LIMIT];
static int16_t cm_modelVertexRemap[CODUO_CM_MODEL_VERTEX_REMAP_LIMIT];

qboolean CM_ComparePoints(const vec3_t left, const vec3_t right);

int32_t CM_LittleShort(int16_t value)
{
    return value;
}

uint32_t CM_LittleLong(uint32_t value)
{
    return value;
}

long double CM_LittleFloat(float value)
{
    return value;
}

uint32_t CM_SignbitsForNormal(const vec3_t normal)
{
    uint32_t signbits = 0;

    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        if (normal[axis] < 0.0f) {
            signbits |= (uint32_t)(1U << axis);
        }
    }

    return signbits;
}

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

    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
#if EMULATE_X87
        /* x87 (DLL 0x0804c019): each edge component is an 80-bit subtract
         * rounded to the float slot (fld;fsub;fstp). */
        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])));
#else
        edge0[axis] = point1[axis] - point0[axis];
        edge1[axis] = point2[axis] - point0[axis];
#endif
    }

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

#if EMULATE_X87
    /* plane[3] = point0 . plane: three 80-bit products and two 80-bit adds,
     * rounded to float once at the store (fmul/faddp x2; fstp). */
    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]))));
#else
    plane[3] = (point0[0] * plane[0]) + (point0[1] * plane[1]) +
               (point0[2] * plane[2]);
#endif
    return qtrue;
}

qboolean CM_NeedsSubdivision(const vec3_t point0, const vec3_t point1,
                             const vec3_t point2, int32_t maxError)
{
    vec3_t bend;

#if EMULATE_X87
    /* x87 (DLL 0x0804c103): bend[axis] = (point0+point2) - (point1+point1) in
     * one 80-bit chain rounded to the float slot (the DLL forms point1+point1
     * with fadd st,st(0)); bendLength = sqrt of the 80-bit dot rounded to a
     * float slot then stored as double for the CRT sqrt (bit-portable), result
     * to a float; the test is maxError < bendLength*0.25f where maxError is
     * loaded exactly by fild and the product is rounded to a float before the
     * compare (fld 0.25f; fmulp; fstp; fild; fld; fucompp). 0.25 is exact, but
     * load it as f32 to mirror fld DWORD. */
    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        x87f p1 = x87f_load_f32(point1[axis]);
        bend[axis] = x87f_store_f32(
            x87f_sub(x87f_add(x87f_load_f32(point0[axis]),
                              x87f_load_f32(point2[axis])),
                     x87f_add(p1, p1)));
    }

    /* The dot is stored DIRECTLY to a double for the CRT sqrt (fstp QWORD),
     * with no intermediate float rounding — unlike VectorNormalize, which
     * rounds lengthSquared to a float slot first. */
    x87f dot = x87f_add(
        x87f_add(x87f_mul(x87f_load_f32(bend[0]), x87f_load_f32(bend[0])),
                 x87f_mul(x87f_load_f32(bend[1]), x87f_load_f32(bend[1]))),
        x87f_mul(x87f_load_f32(bend[2]), x87f_load_f32(bend[2])));
    float bendLength = (float)sqrt(x87f_store_f64(dot));

    float scaled = x87f_store_f32(
        x87f_mul(x87f_load_f32(bendLength), x87f_load_f32(0.25f)));
    return x87f_lt(x87f_load_i32(maxError), x87f_load_f32(scaled))
               ? qtrue
               : qfalse;
#else
    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        bend[axis] = (point0[axis] + point2[axis]) -
                     (point1[axis] + point1[axis]);
    }

    float bendLength =
        (float)sqrt((double)((bend[0] * bend[0]) + (bend[1] * bend[1]) +
                             (bend[2] * bend[2])));

    /* 0x0804c180..0x0804c18a stores the scaled bend to a float slot before
     * comparing it with the converted integer threshold. */
    float scaled = bendLength * 0.25f;
    return (long double)maxError < (long double)scaled ? qtrue : qfalse;
#endif
}

void CM_Subdivide(const vec3_t point0, const vec3_t point1,
                  const vec3_t point2, vec3_t mid01, vec3_t center,
                  vec3_t mid12)
{
#if EMULATE_X87
    /* x87 (DLL 0x0804c1a3): each midpoint is (a+b)*0.5f in one 80-bit chain
     * (fld;fadd;fld 0.5f;fmulp) rounded to the float slot; center reads the
     * just-stored mid01/mid12 floats. 0.5 is exact; loaded as f32 to mirror
     * fld DWORD. */
    const x87f half = x87f_load_f32(0.5f);
    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        mid01[axis] = x87f_store_f32(x87f_mul(
            x87f_add(x87f_load_f32(point0[axis]), x87f_load_f32(point1[axis])),
            half));
        mid12[axis] = x87f_store_f32(x87f_mul(
            x87f_add(x87f_load_f32(point1[axis]), x87f_load_f32(point2[axis])),
            half));
        center[axis] = x87f_store_f32(x87f_mul(
            x87f_add(x87f_load_f32(mid01[axis]), x87f_load_f32(mid12[axis])),
            half));
    }
#else
    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        mid01[axis] = (point0[axis] + point1[axis]) * 0.5f;
        mid12[axis] = (point1[axis] + point2[axis]) * 0.5f;
        center[axis] = (mid01[axis] + mid12[axis]) * 0.5f;
    }
#endif
}

void CM_TransposeGrid(cGrid_t *grid)
{
    if (grid->height < grid->width) {
        for (int32_t x = 0; x < grid->height; ++x) {
            for (int32_t y = x + 1; y < grid->width; ++y) {
                if (y < grid->height) {
                    vec3_t temp = {
                        grid->points[x][y][0],
                        grid->points[x][y][1],
                        grid->points[x][y][2],
                    };

                    grid->points[x][y][0] = grid->points[y][x][0];
                    grid->points[x][y][1] = grid->points[y][x][1];
                    grid->points[x][y][2] = grid->points[y][x][2];
                    grid->points[y][x][0] = temp[0];
                    grid->points[y][x][1] = temp[1];
                    grid->points[y][x][2] = temp[2];
                } else {
                    grid->points[x][y][0] = grid->points[y][x][0];
                    grid->points[x][y][1] = grid->points[y][x][1];
                    grid->points[x][y][2] = grid->points[y][x][2];
                }
            }
        }
    } else {
        for (int32_t x = 0; x < grid->width; ++x) {
            for (int32_t y = x + 1; y < grid->height; ++y) {
                if (y < grid->width) {
                    vec3_t temp = {
                        grid->points[y][x][0],
                        grid->points[y][x][1],
                        grid->points[y][x][2],
                    };

                    grid->points[y][x][0] = grid->points[x][y][0];
                    grid->points[y][x][1] = grid->points[x][y][1];
                    grid->points[y][x][2] = grid->points[x][y][2];
                    grid->points[x][y][0] = temp[0];
                    grid->points[x][y][1] = temp[1];
                    grid->points[x][y][2] = temp[2];
                } else {
                    grid->points[y][x][0] = grid->points[x][y][0];
                    grid->points[y][x][1] = grid->points[x][y][1];
                    grid->points[y][x][2] = grid->points[x][y][2];
                }
            }
        }
    }

    int32_t width = grid->width;
    grid->width = grid->height;
    grid->height = width;

    qboolean wrapWidth = grid->wrapWidth;
    grid->wrapWidth = grid->wrapHeight;
    grid->wrapHeight = wrapWidth;
}

void CM_SetGridWrapWidth(cGrid_t *grid)
{
    int32_t y;

    for (y = 0; y < grid->height; ++y) {
        int32_t axis = 0;
        while (axis < CM_PATCH_VECTOR_AXIS_COUNT) {
#if EMULATE_X87
            /* x87 (DLL 0x0804c8af): delta = points[0] - points[width-1] rounded
             * to the float slot, range-compared against the double wrap epsilon
             * (0.1). */
            float delta = x87f_store_f32(x87f_sub(
                x87f_load_f32(grid->points[0][y][axis]),
                x87f_load_f32(grid->points[grid->width - 1][y][axis])));
            if (x87f_lt(x87f_load_f32(delta),
                        x87f_load_f64(-CODUO_CM_PATCH_WRAP_EPSILON)) ||
                x87f_lt(x87f_load_f64(CODUO_CM_PATCH_WRAP_EPSILON),
                        x87f_load_f32(delta))) {
                break;
            }
#else
            float delta =
                grid->points[0][y][axis] -
                grid->points[grid->width - 1][y][axis];
            if (delta < -CODUO_CM_PATCH_WRAP_EPSILON ||
                CODUO_CM_PATCH_WRAP_EPSILON < delta) {
                break;
            }
#endif
            axis++;
        }

        if (axis != CM_PATCH_VECTOR_AXIS_COUNT) {
            break;
        }
    }

    grid->wrapWidth = y == grid->height ? qtrue : qfalse;
}

void CM_SubdivideGridColumns(cGrid_t *grid, int32_t maxError)
{
    int32_t x = 0;

    while (x < grid->width - 2) {
        int32_t y = 0;
        while (y < grid->height &&
               CM_NeedsSubdivision(grid->points[x][y], grid->points[x + 1][y],
                            grid->points[x + 2][y], maxError) == qfalse) {
            y++;
        }

        if (y == grid->height) {
            x += 2;
            continue;
        }

        /*
         * ORIGINAL_BINARY_BUG: stock shifts the current columns two slots to
         * the right without checking the fixed 129-column point bank. See
         * original-bugs/coduomp-collision-patch-subdivision-grid-overflow.md.
         */
        for (y = 0; y < grid->height; ++y) {
            vec3_t point0 = {
                grid->points[x][y][0],
                grid->points[x][y][1],
                grid->points[x][y][2],
            };
            vec3_t point1 = {
                grid->points[x + 1][y][0],
                grid->points[x + 1][y][1],
                grid->points[x + 1][y][2],
            };
            vec3_t point2 = {
                grid->points[x + 2][y][0],
                grid->points[x + 2][y][1],
                grid->points[x + 2][y][2],
            };

            int32_t source = grid->width;
            while (--source, x + 1 < source) {
                grid->points[source + 2][y][0] = grid->points[source][y][0];
                grid->points[source + 2][y][1] = grid->points[source][y][1];
                grid->points[source + 2][y][2] = grid->points[source][y][2];
            }

            CM_Subdivide(point0, point1, point2, grid->points[x + 1][y],
                         grid->points[x + 2][y], grid->points[x + 3][y]);
        }

        grid->width += 2;
    }
}

void CM_RemoveDegenerateColumns(cGrid_t *grid)
{
    for (int32_t x = 0; x < grid->width - 1; ++x) {
        int32_t y = 0;
        while (y < grid->height &&
               CM_ComparePoints(grid->points[x][y], grid->points[x + 1][y]) !=
                   qfalse) {
            y++;
        }

        if (y != grid->height) {
            continue;
        }

        for (y = 0; y < grid->height; ++y) {
            for (int32_t source = x + 2; source < grid->width; ++source) {
                grid->points[source - 1][y][0] = grid->points[source][y][0];
                grid->points[source - 1][y][1] = grid->points[source][y][1];
                grid->points[source - 1][y][2] = grid->points[source][y][2];
            }
        }

        grid->width--;
        x--;
    }
}

qboolean CM_ComparePoints(const vec3_t left, const vec3_t right)
{
    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
#if EMULATE_X87
        /* x87 (DLL 0x0804cd9c): delta = left - right rounded to the float slot,
         * range-compared against the double point-compare epsilon (0.1). */
        float delta = x87f_store_f32(
            x87f_sub(x87f_load_f32(left[axis]), x87f_load_f32(right[axis])));
        if (x87f_lt(x87f_load_f32(delta),
                    x87f_load_f64(-CODUO_CM_PATCH_POINT_COMPARE_EPSILON)) ||
            x87f_lt(x87f_load_f64(CODUO_CM_PATCH_POINT_COMPARE_EPSILON),
                    x87f_load_f32(delta))) {
            return qfalse;
        }
#else
        float delta = left[axis] - right[axis];
        if (delta < -CODUO_CM_PATCH_POINT_COMPARE_EPSILON ||
            CODUO_CM_PATCH_POINT_COMPARE_EPSILON < delta) {
            return qfalse;
        }
#endif
    }

    return qtrue;
}

qboolean CM_PlaneEqual(const vec4_t left, const vec4_t right,
                       qboolean *flipped)
{
#if EMULATE_X87
    /* x87 (DLL 0x0804d04c): each test is fld left[k]; fsub right[k]; fabs;
     * fld QWORD epsilon; fucompp — i.e. the difference is formed and abs'd in
     * 80-bit and compared against the epsilon widened from a double constant.
     * The `long double` in the source is x87-80-bit on x86 but 128-bit quad on
     * arm64, so route it through the shim to keep the width faithful. The '+'
     * block negates each right[k] via an integer sign-flip (exact for a float)
     * and subtracts, which equals left[k] + right[k] in 80-bit. */
    const x87f normalEps = x87f_load_f64(CODUO_CM_PATCH_PLANE_NORMAL_EPSILON);
    const x87f distEps = x87f_load_f64(CODUO_CM_PATCH_PLANE_DIST_EPSILON);

    if (x87f_lt(x87f_abs(x87f_sub(x87f_load_f32(left[0]),
                                  x87f_load_f32(right[0]))), normalEps) &&
        x87f_lt(x87f_abs(x87f_sub(x87f_load_f32(left[1]),
                                  x87f_load_f32(right[1]))), normalEps) &&
        x87f_lt(x87f_abs(x87f_sub(x87f_load_f32(left[2]),
                                  x87f_load_f32(right[2]))), normalEps) &&
        x87f_lt(x87f_abs(x87f_sub(x87f_load_f32(left[3]),
                                  x87f_load_f32(right[3]))), distEps)) {
        *flipped = qfalse;
        return qtrue;
    }

    if (x87f_lt(x87f_abs(x87f_add(x87f_load_f32(left[0]),
                                  x87f_load_f32(right[0]))), normalEps) &&
        x87f_lt(x87f_abs(x87f_add(x87f_load_f32(left[1]),
                                  x87f_load_f32(right[1]))), normalEps) &&
        x87f_lt(x87f_abs(x87f_add(x87f_load_f32(left[2]),
                                  x87f_load_f32(right[2]))), normalEps) &&
        x87f_lt(x87f_abs(x87f_add(x87f_load_f32(left[3]),
                                  x87f_load_f32(right[3]))), distEps)) {
        *flipped = qtrue;
        return qtrue;
    }

    return qfalse;
#else
    if (fabsl((long double)left[0] - (long double)right[0]) <
            CODUO_CM_PATCH_PLANE_NORMAL_EPSILON &&
        fabsl((long double)left[1] - (long double)right[1]) <
            CODUO_CM_PATCH_PLANE_NORMAL_EPSILON &&
        fabsl((long double)left[2] - (long double)right[2]) <
            CODUO_CM_PATCH_PLANE_NORMAL_EPSILON &&
        fabsl((long double)left[3] - (long double)right[3]) <
            CODUO_CM_PATCH_PLANE_DIST_EPSILON) {
        *flipped = qfalse;
        return qtrue;
    }

    if (fabsl((long double)left[0] + (long double)right[0]) <
            CODUO_CM_PATCH_PLANE_NORMAL_EPSILON &&
        fabsl((long double)left[1] + (long double)right[1]) <
            CODUO_CM_PATCH_PLANE_NORMAL_EPSILON &&
        fabsl((long double)left[2] + (long double)right[2]) <
            CODUO_CM_PATCH_PLANE_NORMAL_EPSILON &&
        fabsl((long double)left[3] + (long double)right[3]) <
            CODUO_CM_PATCH_PLANE_DIST_EPSILON) {
        *flipped = qtrue;
        return qtrue;
    }

    return qfalse;
#endif
}

void CM_SnapVector(vec3_t normal)
{
#if EMULATE_X87
    /* x87 (DLL 0x0804d1a1): |normal[axis] - 1| via fld;fld1;fsubp;fabs and
     * |normal[axis] + 1| via fld;fld -1.0f;fsubp;fabs, each compared in 80-bit
     * against the double snap epsilon. 1.0 is exact in every format, so the
     * shim loads it as a double constant. The source's `long double` is not
     * 80-bit on arm64 — route through the shim. */
    const x87f one = x87f_load_f64(1.0);
    const x87f snapEps = x87f_load_f64(CODUO_CM_PATCH_SNAP_VECTOR_EPSILON);
    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        if (x87f_lt(x87f_abs(x87f_sub(x87f_load_f32(normal[axis]), one)),
                    snapEps)) {
            normal[0] = 0.0f;
            normal[1] = 0.0f;
            normal[2] = 0.0f;
            normal[axis] = 1.0f;
            return;
        }

        if (x87f_lt(x87f_abs(x87f_add(x87f_load_f32(normal[axis]), one)),
                    snapEps)) {
            normal[0] = 0.0f;
            normal[1] = 0.0f;
            normal[2] = 0.0f;
            normal[axis] = -1.0f;
            return;
        }
    }
#else
    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        if (fabsl((long double)normal[axis] - 1.0L) <
            CODUO_CM_PATCH_SNAP_VECTOR_EPSILON) {
            normal[0] = 0.0f;
            normal[1] = 0.0f;
            normal[2] = 0.0f;
            normal[axis] = 1.0f;
            return;
        }

        if (fabsl((long double)normal[axis] + 1.0L) <
            CODUO_CM_PATCH_SNAP_VECTOR_EPSILON) {
            normal[0] = 0.0f;
            normal[1] = 0.0f;
            normal[2] = 0.0f;
            normal[axis] = -1.0f;
            return;
        }
    }
#endif
}

int32_t CM_FindPlane(const vec4_t plane, qboolean *flipped)
{
    for (int32_t planeIndex = 0; planeIndex < cm_patchPlaneCount;
         ++planeIndex) {
        vec4_t storedPlane = {
            cm_patchPlanes[planeIndex].normal[0],
            cm_patchPlanes[planeIndex].normal[1],
            cm_patchPlanes[planeIndex].normal[2],
            cm_patchPlanes[planeIndex].dist,
        };

        if (CM_PlaneEqual(storedPlane, plane, flipped) != qfalse) {
            return planeIndex;
        }
    }

    if (cm_patchPlaneCount == CM_PATCH_MAX_PLANES) {
        Com_Error(1, "\x15" "MAX_PATCH_PLANES");
    }

    int32_t planeIndex = cm_patchPlaneCount;
    cm_patchPlanes[planeIndex].normal[0] = plane[0];
    cm_patchPlanes[planeIndex].normal[1] = plane[1];
    cm_patchPlanes[planeIndex].normal[2] = plane[2];
    cm_patchPlanes[planeIndex].dist = plane[3];
    cm_patchPlanes[planeIndex].signbits = CM_SignbitsForNormal(plane);
    cm_patchPlaneCount++;

    *flipped = qfalse;
    return planeIndex;
}

int32_t CM_FindPlane2(const vec3_t point0, const vec3_t point1,
                      const vec3_t point2)
{
    vec4_t plane;

    if (CM_PlaneFromPoints(plane, point0, point1, point2) == qfalse) {
        return -1;
    }

    for (int32_t planeIndex = 0; planeIndex < cm_patchPlaneCount;
         ++planeIndex) {
        const patchPlane_t *storedPlane =
            &cm_patchPlanes[planeIndex];

#if EMULATE_X87
        /* x87 (DLL 0x0804d3c2). Sign test: plane . storedPlane->normal formed
         * entirely in 80-bit and compared against zero (fmul/faddp x2; fldz;
         * fucompp) — NO rounding store, so it stays 80-bit (the source's
         * `long double`, which is 128-bit quad on arm64). */
        x87f signDot = x87f_add(
            x87f_add(x87f_mul(x87f_load_f32(plane[0]),
                              x87f_load_f32(storedPlane->normal[0])),
                     x87f_mul(x87f_load_f32(plane[1]),
                              x87f_load_f32(storedPlane->normal[1]))),
            x87f_mul(x87f_load_f32(plane[2]),
                     x87f_load_f32(storedPlane->normal[2])));
        if (x87f_lt(signDot, x87f_load_f64(0.0))) {
            continue;
        }

        /* Each pointNDist: point . normal - dist in 80-bit, rounded to the
         * float slot (fstp), then the reloaded float is range-compared against
         * the double point-plane epsilon (fld QWORD ±0.1; fucompp). */
        const x87f pointPlaneEps =
            x87f_load_f64(CODUO_CM_PATCH_POINT_PLANE_EPSILON);
        const x87f pointPlaneNegEps =
            x87f_load_f64(-CODUO_CM_PATCH_POINT_PLANE_EPSILON);

        float point0Dist = x87f_store_f32(x87f_sub(
            x87f_add(x87f_add(x87f_mul(x87f_load_f32(point0[0]),
                                       x87f_load_f32(storedPlane->normal[0])),
                              x87f_mul(x87f_load_f32(point0[1]),
                                       x87f_load_f32(storedPlane->normal[1]))),
                     x87f_mul(x87f_load_f32(point0[2]),
                              x87f_load_f32(storedPlane->normal[2]))),
            x87f_load_f32(storedPlane->dist)));
        if (x87f_lt(x87f_load_f32(point0Dist), pointPlaneNegEps) ||
            x87f_lt(pointPlaneEps, x87f_load_f32(point0Dist))) {
            continue;
        }

        float point1Dist = x87f_store_f32(x87f_sub(
            x87f_add(x87f_add(x87f_mul(x87f_load_f32(point1[0]),
                                       x87f_load_f32(storedPlane->normal[0])),
                              x87f_mul(x87f_load_f32(point1[1]),
                                       x87f_load_f32(storedPlane->normal[1]))),
                     x87f_mul(x87f_load_f32(point1[2]),
                              x87f_load_f32(storedPlane->normal[2]))),
            x87f_load_f32(storedPlane->dist)));
        if (x87f_lt(x87f_load_f32(point1Dist), pointPlaneNegEps) ||
            x87f_lt(pointPlaneEps, x87f_load_f32(point1Dist))) {
            continue;
        }

        float point2Dist = x87f_store_f32(x87f_sub(
            x87f_add(x87f_add(x87f_mul(x87f_load_f32(point2[0]),
                                       x87f_load_f32(storedPlane->normal[0])),
                              x87f_mul(x87f_load_f32(point2[1]),
                                       x87f_load_f32(storedPlane->normal[1]))),
                     x87f_mul(x87f_load_f32(point2[2]),
                              x87f_load_f32(storedPlane->normal[2]))),
            x87f_load_f32(storedPlane->dist)));
        if (x87f_lt(x87f_load_f32(point2Dist), pointPlaneNegEps) ||
            x87f_lt(pointPlaneEps, x87f_load_f32(point2Dist))) {
            continue;
        }

        return planeIndex;
#else
        if ((long double)plane[0] * (long double)storedPlane->normal[0] +
                (long double)plane[1] * (long double)storedPlane->normal[1] +
                (long double)plane[2] * (long double)storedPlane->normal[2] <
            0.0) {
            continue;
        }

        float point0Dist =
            (point0[0] * storedPlane->normal[0]) +
            (point0[1] * storedPlane->normal[1]) +
            (point0[2] * storedPlane->normal[2]) - storedPlane->dist;
        if (point0Dist < -CODUO_CM_PATCH_POINT_PLANE_EPSILON ||
            CODUO_CM_PATCH_POINT_PLANE_EPSILON < point0Dist) {
            continue;
        }

        float point1Dist =
            (point1[0] * storedPlane->normal[0]) +
            (point1[1] * storedPlane->normal[1]) +
            (point1[2] * storedPlane->normal[2]) - storedPlane->dist;
        if (point1Dist < -CODUO_CM_PATCH_POINT_PLANE_EPSILON ||
            CODUO_CM_PATCH_POINT_PLANE_EPSILON < point1Dist) {
            continue;
        }

        float point2Dist =
            (point2[0] * storedPlane->normal[0]) +
            (point2[1] * storedPlane->normal[1]) +
            (point2[2] * storedPlane->normal[2]) - storedPlane->dist;

        if (point2Dist < -CODUO_CM_PATCH_POINT_PLANE_EPSILON ||
            CODUO_CM_PATCH_POINT_PLANE_EPSILON < point2Dist) {
            continue;
        }

        return planeIndex;
#endif
    }

    if (cm_patchPlaneCount == CM_PATCH_MAX_PLANES) {
        Com_Error(1, "\x15" "MAX_PATCH_PLANES");
    }

    int32_t planeIndex = cm_patchPlaneCount;
    cm_patchPlanes[planeIndex].normal[0] = plane[0];
    cm_patchPlanes[planeIndex].normal[1] = plane[1];
    cm_patchPlanes[planeIndex].normal[2] = plane[2];
    cm_patchPlanes[planeIndex].dist = plane[3];
    cm_patchPlanes[planeIndex].signbits = CM_SignbitsForNormal(plane);
    cm_patchPlaneCount++;

    return planeIndex;
}

int32_t CM_PointOnPlaneSide(const vec3_t point, int32_t planeIndex)
{
    if (planeIndex == -1) {
        return CODUO_CM_PATCH_PLANE_SIDE_ON;
    }

    const patchPlane_t *plane = &cm_patchPlanes[planeIndex];
#if EMULATE_X87
    /* x87 (DLL 0x0804d6f7): point . normal - dist in one 80-bit chain rounded
     * to the float dist slot, then range-compared against the double
     * point-plane epsilon (0.1). */
    float dist = x87f_store_f32(x87f_sub(
        x87f_add(x87f_add(x87f_mul(x87f_load_f32(point[0]),
                                   x87f_load_f32(plane->normal[0])),
                          x87f_mul(x87f_load_f32(point[1]),
                                   x87f_load_f32(plane->normal[1]))),
                 x87f_mul(x87f_load_f32(point[2]),
                          x87f_load_f32(plane->normal[2]))),
        x87f_load_f32(plane->dist)));

    if (x87f_lt(x87f_load_f64(CODUO_CM_PATCH_POINT_PLANE_EPSILON),
                x87f_load_f32(dist))) {
        return CODUO_CM_PATCH_PLANE_SIDE_FRONT;
    }

    if (x87f_lt(x87f_load_f32(dist),
                x87f_load_f64(-CODUO_CM_PATCH_POINT_PLANE_EPSILON))) {
        return CODUO_CM_PATCH_PLANE_SIDE_BACK;
    }
#else
    float dist = (point[0] * plane->normal[0]) +
                 (point[1] * plane->normal[1]) +
                 (point[2] * plane->normal[2]) - plane->dist;

    if (CODUO_CM_PATCH_POINT_PLANE_EPSILON < dist) {
        return CODUO_CM_PATCH_PLANE_SIDE_FRONT;
    }

    if (dist < -CODUO_CM_PATCH_POINT_PLANE_EPSILON) {
        return CODUO_CM_PATCH_PLANE_SIDE_BACK;
    }
#endif

    return CODUO_CM_PATCH_PLANE_SIDE_ON;
}

int32_t CM_GridPlane(
    int32_t planeGrid[CODUO_CM_PATCH_PLANE_GRID_SIZE]
                     [CODUO_CM_PATCH_PLANE_GRID_SIZE][2],
    int32_t row, int32_t column, int32_t planeSlot)
{
    int32_t planeIndex = planeGrid[row][column][planeSlot];
    if (planeIndex == -1) {
        planeIndex = planeGrid[row][column][planeSlot == 0 ? 1 : 0];
    }

    return planeIndex;
}

int32_t CM_EdgePlaneNum(
    const cGrid_t *grid,
    int32_t planeGrid[CODUO_CM_PATCH_PLANE_GRID_SIZE]
                     [CODUO_CM_PATCH_PLANE_GRID_SIZE][2],
    int32_t x, int32_t y, int32_t edgeName)
{
    const float *point0;
    const float *point1;
    int32_t planeSlot;
    qboolean reversePoints;

    switch (edgeName) {
    case CODUO_CM_PATCH_EDGE_BOTTOM:
        point0 = grid->points[x][y];
        point1 = grid->points[x + 1][y];
        planeSlot = 0;
        reversePoints = qfalse;
        break;
    case CODUO_CM_PATCH_EDGE_RIGHT:
        point0 = grid->points[x + 1][y];
        point1 = grid->points[x + 1][y + 1];
        planeSlot = 0;
        reversePoints = qfalse;
        break;
    case CODUO_CM_PATCH_EDGE_TOP:
        point0 = grid->points[x][y + 1];
        point1 = grid->points[x + 1][y + 1];
        planeSlot = 1;
        reversePoints = qtrue;
        break;
    case CODUO_CM_PATCH_EDGE_LEFT:
        point0 = grid->points[x][y];
        point1 = grid->points[x][y + 1];
        planeSlot = 1;
        reversePoints = qtrue;
        break;
    case CODUO_CM_PATCH_EDGE_DIAGONAL_DESCENDING:
        point0 = grid->points[x + 1][y + 1];
        point1 = grid->points[x][y];
        planeSlot = 0;
        reversePoints = qfalse;
        break;
    case CODUO_CM_PATCH_EDGE_DIAGONAL_ASCENDING:
        point0 = grid->points[x][y];
        point1 = grid->points[x + 1][y + 1];
        planeSlot = 1;
        reversePoints = qfalse;
        break;
    default:
        Com_Error(1, "\x15" "CM_EdgePlaneNum: bad edge name");
        return -1;
    }

    int32_t adjacentPlane = CM_GridPlane(planeGrid, x, y, planeSlot);
    if (adjacentPlane < 0) {
        return -1;
    }

#if EMULATE_X87
    /* x87 (DLL 0x0804d8b0): offsetPoint = point0 + normal * PUSH_DISTANCE,
     * each axis one 80-bit chain (fld normal; fld 4.0f; fmulp; fld point0;
     * faddp) rounded to the float slot. PUSH_DISTANCE=4 loads as float 4.0. */
    vec3_t offsetPoint;
    for (int32_t k = 0; k < CM_PATCH_VECTOR_AXIS_COUNT; ++k) {
        offsetPoint[k] = x87f_store_f32(x87f_add(
            x87f_load_f32(point0[k]),
            x87f_mul(x87f_load_f32(cm_patchPlanes[adjacentPlane].normal[k]),
                     x87f_load_f32(
                         (float)CODUO_CM_PATCH_EDGE_PLANE_PUSH_DISTANCE))));
    }
#else
    vec3_t offsetPoint = {
        point0[0] + cm_patchPlanes[adjacentPlane].normal[0] *
                        CODUO_CM_PATCH_EDGE_PLANE_PUSH_DISTANCE,
        point0[1] + cm_patchPlanes[adjacentPlane].normal[1] *
                        CODUO_CM_PATCH_EDGE_PLANE_PUSH_DISTANCE,
        point0[2] + cm_patchPlanes[adjacentPlane].normal[2] *
                        CODUO_CM_PATCH_EDGE_PLANE_PUSH_DISTANCE,
    };
#endif

    if (reversePoints != qfalse) {
        return CM_FindPlane2(point1, point0, offsetPoint);
    }

    return CM_FindPlane2(point0, point1, offsetPoint);
}

void CM_SetBorderInward(
    facet_t *facet, const cGrid_t *grid,
    int32_t planeGrid[CODUO_CM_PATCH_PLANE_GRID_SIZE]
                     [CODUO_CM_PATCH_PLANE_GRID_SIZE][2],
    int32_t x, int32_t y, int32_t triangleMode)
{
    (void)planeGrid;

    const float *vertices[4];
    int32_t vertexCount;

    if (triangleMode == CODUO_CM_PATCH_BORDER_INWARD_LOWER_TRIANGLE) {
        vertices[0] = grid->points[x][y];
        vertices[1] = grid->points[x + 1][y];
        vertices[2] = grid->points[x + 1][y + 1];
        vertexCount = 3;
    } else if (triangleMode == CODUO_CM_PATCH_BORDER_INWARD_QUAD) {
        vertices[0] = grid->points[x][y];
        vertices[1] = grid->points[x + 1][y];
        vertices[2] = grid->points[x + 1][y + 1];
        vertices[3] = grid->points[x][y + 1];
        vertexCount = 4;
    } else if (triangleMode == CODUO_CM_PATCH_BORDER_INWARD_UPPER_TRIANGLE) {
        vertices[0] = grid->points[x + 1][y + 1];
        vertices[1] = grid->points[x][y + 1];
        vertices[2] = grid->points[x][y];
        vertexCount = 3;
    } else {
        Com_Error(0, "\x15" "CM_SetBorderInward: bad parameter");
        vertexCount = 0;
    }

    for (int32_t planeIndex = 0; planeIndex < facet->numBorders;
         ++planeIndex) {
        int32_t frontCount = 0;
        int32_t backCount = 0;

        for (int32_t vertexIndex = 0; vertexIndex < vertexCount;
             ++vertexIndex) {
            int32_t side = CM_PointOnPlaneSide(vertices[vertexIndex],
                                        facet->borderPlanes[planeIndex]);
            if (side == CODUO_CM_PATCH_PLANE_SIDE_FRONT) {
                frontCount++;
            } else if (side == CODUO_CM_PATCH_PLANE_SIDE_BACK) {
                backCount++;
            }
        }

        if (frontCount != 0 && backCount == 0) {
            facet->borderInward[planeIndex] = qtrue;
        } else if (backCount != 0 && frontCount == 0) {
            facet->borderInward[planeIndex] = qfalse;
        } else if (frontCount == 0 && backCount == 0) {
            facet->borderPlanes[planeIndex] = -1;
        } else {
            Com_DPrintf(
                "WARNING: CM_SetBorderInward: mixed plane sides\n");
            facet->borderInward[planeIndex] = qfalse;
        }
    }
}

void PrintWinding(const winding_t *winding)
{
    for (int32_t pointIndex = 0; pointIndex < winding->numpoints;
         ++pointIndex) {
        printf("(%5.1f, %5.1f, %5.1f)\n",
               (double)winding->p[pointIndex][0],
               (double)winding->p[pointIndex][1],
               (double)winding->p[pointIndex][2]);
    }
}

winding_t *AllocWinding(int32_t pointCapacity)
{
    cm_windingActiveCount++;
    if (cm_windingPeakActiveCount < cm_windingActiveCount) {
        cm_windingPeakActiveCount = cm_windingActiveCount;
    }

    /* 0x080514a7..0x080514b6 forms 4 + pointCapacity*12 in a dword. */
    uint32_t windingBytes =
        (uint32_t)sizeof(winding_t) +
        (uint32_t)pointCapacity *
            (uint32_t)sizeof(((winding_t *)0)->p[0]);
    return Z_Malloc((size_t)windingBytes);
}

void FreeWinding(winding_t *winding)
{
    if (winding->numpoints == CODUO_CM_WINDING_FREED_SENTINEL) {
        Com_Error(0, "\x15" "FreeWinding: freed a freed winding");
    }

    winding->numpoints = CODUO_CM_WINDING_FREED_SENTINEL;
    cm_windingActiveCount--;
    Z_Free(winding);
}

void RemoveColinearPoints(winding_t *winding)
{
    vec3_t keptPoints[CODUO_CM_WINDING_TEMP_POINT_LIMIT];
    int32_t keptCount = 0;

    for (int32_t pointIndex = 0; pointIndex < winding->numpoints;
         ++pointIndex) {
        int32_t nextIndex = (pointIndex + 1) % winding->numpoints;
        int32_t previousIndex =
            (winding->numpoints + pointIndex - 1) % winding->numpoints;
        vec3_t nextEdge = {
            winding->p[nextIndex][0] - winding->p[pointIndex][0],
            winding->p[nextIndex][1] - winding->p[pointIndex][1],
            winding->p[nextIndex][2] - winding->p[pointIndex][2],
        };
        vec3_t previousEdge = {
            winding->p[pointIndex][0] -
                winding->p[previousIndex][0],
            winding->p[pointIndex][1] -
                winding->p[previousIndex][1],
            winding->p[pointIndex][2] -
                winding->p[previousIndex][2],
        };

        VectorNormalize2(nextEdge, nextEdge);
        VectorNormalize2(previousEdge, previousEdge);

        long double dot = ((long double)nextEdge[0] * previousEdge[0]) +
                          ((long double)nextEdge[1] * previousEdge[1]) +
                          ((long double)nextEdge[2] * previousEdge[2]);
        if (dot < CODUO_CM_WINDING_COLINEAR_DOT_EPSILON) {
            keptPoints[keptCount][0] = winding->p[pointIndex][0];
            keptPoints[keptCount][1] = winding->p[pointIndex][1];
            keptPoints[keptCount][2] = winding->p[pointIndex][2];
            keptCount++;
        }
    }

    if (keptCount != winding->numpoints) {
        winding->numpoints = keptCount;
        uint32_t copyBytes =
            (uint32_t)keptCount * (uint32_t)sizeof(keptPoints[0]);
        memcpy(winding->p, keptPoints, (size_t)copyBytes);
    }
}

void WindingPlane(const winding_t *winding, vec3_t normal,
                  float *dist)
{
    vec3_t edge0 = {
        winding->p[1][0] - winding->p[0][0],
        winding->p[1][1] - winding->p[0][1],
        winding->p[1][2] - winding->p[0][2],
    };
    vec3_t edge1 = {
        winding->p[2][0] - winding->p[0][0],
        winding->p[2][1] - winding->p[0][1],
        winding->p[2][2] - winding->p[0][2],
    };

    CrossProduct(edge1, edge0, normal);
    VectorNormalize2(normal, normal);
    *dist = (winding->p[0][0] * normal[0]) +
            (winding->p[0][1] * normal[1]) +
            (winding->p[0][2] * normal[2]);
}

float WindingArea(const winding_t *winding)
{
    float area = 0.0f;

    for (int32_t pointIndex = 2; pointIndex < winding->numpoints;
         ++pointIndex) {
        vec3_t edge0 = {
            winding->p[pointIndex - 1][0] - winding->p[0][0],
            winding->p[pointIndex - 1][1] - winding->p[0][1],
            winding->p[pointIndex - 1][2] - winding->p[0][2],
        };
        vec3_t edge1 = {
            winding->p[pointIndex][0] - winding->p[0][0],
            winding->p[pointIndex][1] - winding->p[0][1],
            winding->p[pointIndex][2] - winding->p[0][2],
        };
        vec3_t cross;

        CrossProduct(edge0, edge1, cross);
        area +=
            (float)sqrt((double)((cross[0] * cross[0]) +
                                 (cross[1] * cross[1]) +
                                 (cross[2] * cross[2]))) *
            0.5f;
    }

    return area;
}

void WindingBounds(const winding_t *winding, vec3_t mins, vec3_t maxs)
{
    mins[2] = (float)CODUO_CM_WINDING_WORLD_COORD_LIMIT;
    mins[1] = (float)CODUO_CM_WINDING_WORLD_COORD_LIMIT;
    mins[0] = (float)CODUO_CM_WINDING_WORLD_COORD_LIMIT;
    maxs[2] = (float)-CODUO_CM_WINDING_WORLD_COORD_LIMIT;
    maxs[1] = (float)-CODUO_CM_WINDING_WORLD_COORD_LIMIT;
    maxs[0] = (float)-CODUO_CM_WINDING_WORLD_COORD_LIMIT;

    for (int32_t pointIndex = 0; pointIndex < winding->numpoints;
         ++pointIndex) {
        for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT;
             ++axis) {
            if (winding->p[pointIndex][axis] < mins[axis]) {
                mins[axis] = winding->p[pointIndex][axis];
            }
            if (maxs[axis] < winding->p[pointIndex][axis]) {
                maxs[axis] = winding->p[pointIndex][axis];
            }
        }
    }
}

void WindingCenter(const winding_t *winding, vec3_t center)
{
    center[0] = 0.0f;
    center[1] = 0.0f;
    center[2] = 0.0f;

    for (int32_t pointIndex = 0; pointIndex < winding->numpoints;
         ++pointIndex) {
        center[0] += winding->p[pointIndex][0];
        center[1] += winding->p[pointIndex][1];
        center[2] += winding->p[pointIndex][2];
    }

    /* 0x8051b45: bare fild feeds the divide; no (float) cast rounding. */
    float scale = 1.0f / winding->numpoints;
    center[0] *= scale;
    center[1] *= scale;
    center[2] *= scale;
}

winding_t *BaseWindingForPlane(const vec3_t normal, float dist)
{
    int32_t majorAxis = -1;
    float majorAxisMagnitude = (float)-CODUO_CM_WINDING_WORLD_COORD_LIMIT;

    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        float magnitude = fabsf(normal[axis]);
        if (majorAxisMagnitude < magnitude) {
            majorAxis = axis;
            majorAxisMagnitude = magnitude;
        }
    }

    if (majorAxis == -1) {
        Com_Error(1, "\x15" "BaseWindingForPlane: no axis found");
    }

    vec3_t up = {0.0f, 0.0f, 0.0f};
    if (majorAxis < 2) {
        up[2] = 1.0f;
    } else {
        up[0] = 1.0f;
    }

#if EMULATE_X87
    /* x87 (DLL 0x08051b87). projection = up . normal (80-bit dot, one store);
     * up -= projection*normal (fld projection; fchs; fmul normal; fld up;
     * faddp; store — i.e. up + (-projection)*normal in 80-bit); origin =
     * normal*dist; up,right *= LIMIT. VectorNormalize2/CrossProduct are the
     * emulated variants. */
    float projection = x87f_store_f32(x87f_add(
        x87f_add(x87f_mul(x87f_load_f32(up[0]), x87f_load_f32(normal[0])),
                 x87f_mul(x87f_load_f32(up[1]), x87f_load_f32(normal[1]))),
        x87f_mul(x87f_load_f32(up[2]), x87f_load_f32(normal[2]))));
    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        up[axis] = x87f_store_f32(
            x87f_add(x87f_load_f32(up[axis]),
                     x87f_mul(x87f_neg(x87f_load_f32(projection)),
                              x87f_load_f32(normal[axis]))));
    }
    VectorNormalize2(up, up);

    vec3_t origin;
    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        origin[axis] = x87f_store_f32(
            x87f_mul(x87f_load_f32(normal[axis]), x87f_load_f32(dist)));
    }
    vec3_t right;
    CrossProduct(up, normal, right);

    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        up[axis] = x87f_store_f32(x87f_mul(
            x87f_load_f32(up[axis]),
            x87f_load_f32((float)CODUO_CM_WINDING_WORLD_COORD_LIMIT)));
        right[axis] = x87f_store_f32(x87f_mul(
            x87f_load_f32(right[axis]),
            x87f_load_f32((float)CODUO_CM_WINDING_WORLD_COORD_LIMIT)));
    }

    winding_t *winding = AllocWinding(4);
    /* Each point is TWO 80-bit chains with the first rounded to a float slot
     * (the DLL stores origin±right to the point slot, reloads it, then ±up):
     * points = fstps(fstps(origin±right) ±up). */
    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        float originMinusRight = x87f_store_f32(
            x87f_sub(x87f_load_f32(origin[axis]), x87f_load_f32(right[axis])));
        float originPlusRight = x87f_store_f32(
            x87f_add(x87f_load_f32(origin[axis]), x87f_load_f32(right[axis])));
        winding->p[0][axis] = x87f_store_f32(x87f_add(
            x87f_load_f32(originMinusRight), x87f_load_f32(up[axis])));
        winding->p[1][axis] = x87f_store_f32(x87f_add(
            x87f_load_f32(originPlusRight), x87f_load_f32(up[axis])));
        winding->p[2][axis] = x87f_store_f32(x87f_sub(
            x87f_load_f32(originPlusRight), x87f_load_f32(up[axis])));
        winding->p[3][axis] = x87f_store_f32(x87f_sub(
            x87f_load_f32(originMinusRight), x87f_load_f32(up[axis])));
    }
    winding->numpoints = 4;

    return winding;
#else
    float projection =
        (up[0] * normal[0]) + (up[1] * normal[1]) + (up[2] * normal[2]);
    up[0] -= projection * normal[0];
    up[1] -= projection * normal[1];
    up[2] -= projection * normal[2];
    VectorNormalize2(up, up);

    vec3_t origin = {
        normal[0] * dist,
        normal[1] * dist,
        normal[2] * dist,
    };
    vec3_t right;
    CrossProduct(up, normal, right);

    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        up[axis] *= (float)CODUO_CM_WINDING_WORLD_COORD_LIMIT;
        right[axis] *= (float)CODUO_CM_WINDING_WORLD_COORD_LIMIT;
    }

    winding_t *winding = AllocWinding(4);
    /* The original rounds the origin±right subterm to a float before adding
     * ±up (it spills to the point slot and reloads), so use explicit float
     * temporaries: a natural `origin - right + up` would instead keep the
     * subterm in an 80-bit register on both x87 and -mfpmath=387 and not
     * match the stock binary. */
    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        float originMinusRight = origin[axis] - right[axis];
        float originPlusRight = origin[axis] + right[axis];
        winding->p[0][axis] = originMinusRight + up[axis];
        winding->p[1][axis] = originPlusRight + up[axis];
        winding->p[2][axis] = originPlusRight - up[axis];
        winding->p[3][axis] = originMinusRight - up[axis];
    }
    winding->numpoints = 4;

    return winding;
#endif
}

winding_t *CopyWinding(const winding_t *winding)
{
    winding_t *copy = AllocWinding(winding->numpoints);
    /* 0x08051ee3..0x08051efb passes the dword 4 + numpoints*12 to memcpy. */
    uint32_t copyBytes =
        (uint32_t)sizeof(*winding) +
        (uint32_t)winding->numpoints *
            (uint32_t)sizeof(winding->p[0]);
    memcpy(copy, winding, (size_t)copyBytes);
    return copy;
}

winding_t *ReverseWinding(const winding_t *winding)
{
    winding_t *reverse = AllocWinding(winding->numpoints);

    for (int32_t pointIndex = 0; pointIndex < winding->numpoints;
         ++pointIndex) {
        int32_t sourceIndex = winding->numpoints - pointIndex - 1;
        reverse->p[pointIndex][0] = winding->p[sourceIndex][0];
        reverse->p[pointIndex][1] = winding->p[sourceIndex][1];
        reverse->p[pointIndex][2] = winding->p[sourceIndex][2];
    }

    reverse->numpoints = winding->numpoints;
    return reverse;
}

void ClipWindingEpsilon(const winding_t *winding,
                        const vec3_t normal, float dist, float epsilon,
                        winding_t **front,
                        winding_t **back)
{
    int32_t sideCounts[3] = {0, 0, 0};
    int32_t sides[CODUO_CM_WINDING_TEMP_POINT_LIMIT];
    float dists[CODUO_CM_WINDING_TEMP_POINT_LIMIT];

    for (int32_t pointIndex = 0; pointIndex < winding->numpoints;
         ++pointIndex) {
#if EMULATE_X87
        /* classification dot: points . normal in 80-bit rounded to the float
         * global, then a second store for -= dist (matches ChopWindingInPlace
         * and DLL 0x08051fe8). */
        x87f cd = x87f_mul(x87f_load_f32(winding->p[pointIndex][0]),
                           x87f_load_f32(normal[0]));
        cd = x87f_add(cd, x87f_mul(x87f_load_f32(winding->p[pointIndex][1]),
                                   x87f_load_f32(normal[1])));
        cd = x87f_add(cd, x87f_mul(x87f_load_f32(winding->p[pointIndex][2]),
                                   x87f_load_f32(normal[2])));
        cm_windingSplitDist = x87f_store_f32(cd);
        cm_windingSplitDist =
            x87f_store_f32(x87f_sub(x87f_load_f32(cm_windingSplitDist),
                                    x87f_load_f32(dist)));
#else
        cm_windingSplitDist =
            (winding->p[pointIndex][0] * normal[0]) +
            (winding->p[pointIndex][1] * normal[1]) +
            (winding->p[pointIndex][2] * normal[2]);
        cm_windingSplitDist -= dist;
#endif
        dists[pointIndex] = cm_windingSplitDist;

        if (epsilon < cm_windingSplitDist) {
            sides[pointIndex] = CODUO_CM_PATCH_PLANE_SIDE_FRONT;
        } else if (cm_windingSplitDist < -epsilon) {
            sides[pointIndex] = CODUO_CM_PATCH_PLANE_SIDE_BACK;
        } else {
            sides[pointIndex] = CODUO_CM_PATCH_PLANE_SIDE_ON;
        }
        sideCounts[sides[pointIndex]]++;
    }

    sides[winding->numpoints] = sides[0];
    dists[winding->numpoints] = dists[0];
    *back = NULL;
    *front = NULL;

    if (sideCounts[CODUO_CM_PATCH_PLANE_SIDE_FRONT] == 0) {
        *back = CopyWinding(winding);
        return;
    }

    if (sideCounts[CODUO_CM_PATCH_PLANE_SIDE_BACK] == 0) {
        *front = CopyWinding(winding);
        return;
    }

    int32_t splitCapacity =
        winding->numpoints + CODUO_CM_WINDING_SPLIT_EXTRA_POINTS;
    winding_t *frontWinding = AllocWinding(splitCapacity);
    *front = frontWinding;
    winding_t *backWinding = AllocWinding(splitCapacity);
    *back = backWinding;

    for (int32_t pointIndex = 0; pointIndex < winding->numpoints;
         ++pointIndex) {
        const float *point = winding->p[pointIndex];

        if (sides[pointIndex] == CODUO_CM_PATCH_PLANE_SIDE_ON) {
            for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT;
                 ++axis) {
                frontWinding->p[frontWinding->numpoints][axis] =
                    point[axis];
                backWinding->p[backWinding->numpoints][axis] =
                    point[axis];
            }
            frontWinding->numpoints++;
            backWinding->numpoints++;
        } else {
            if (sides[pointIndex] == CODUO_CM_PATCH_PLANE_SIDE_FRONT) {
                for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT;
                     ++axis) {
                    frontWinding->p[frontWinding->numpoints][axis] =
                        point[axis];
                }
                frontWinding->numpoints++;
            }
            if (sides[pointIndex] == CODUO_CM_PATCH_PLANE_SIDE_BACK) {
                for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT;
                     ++axis) {
                    backWinding->p[backWinding->numpoints][axis] =
                        point[axis];
                }
                backWinding->numpoints++;
            }

            if (sides[pointIndex + 1] != CODUO_CM_PATCH_PLANE_SIDE_ON &&
                sides[pointIndex + 1] != sides[pointIndex]) {
                const float *nextPoint =
                    winding->p[(pointIndex + 1) % winding->numpoints];
#if EMULATE_X87
                /* split ratio dists[i]/(dists[i]-dists[i+1]) in 80-bit rounded
                 * to the float global (matches ChopWindingInPlace). */
                cm_windingSplitDist = x87f_store_f32(
                    x87f_div(x87f_load_f32(dists[pointIndex]),
                             x87f_sub(x87f_load_f32(dists[pointIndex]),
                                      x87f_load_f32(dists[pointIndex + 1]))));
#else
                cm_windingSplitDist =
                    dists[pointIndex] /
                    (dists[pointIndex] - dists[pointIndex + 1]);
#endif

                vec3_t splitPoint;
                for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT;
                     ++axis) {
                    if (normal[axis] == CODUO_CM_WINDING_AXIS_NORMAL) {
                        splitPoint[axis] = dist;
                    } else if (normal[axis] == -CODUO_CM_WINDING_AXIS_NORMAL) {
                        splitPoint[axis] = -dist;
                    } else {
#if EMULATE_X87
                        /* lerp point + (nextPoint-point)*ratio in 80-bit
                         * rounded to the float slot. */
                        splitPoint[axis] = x87f_store_f32(x87f_add(
                            x87f_load_f32(point[axis]),
                            x87f_mul(x87f_sub(x87f_load_f32(nextPoint[axis]),
                                              x87f_load_f32(point[axis])),
                                     x87f_load_f32(cm_windingSplitDist))));
#else
                        splitPoint[axis] =
                            point[axis] +
                            (nextPoint[axis] - point[axis]) *
                                cm_windingSplitDist;
#endif
                    }
                }

                for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT;
                     ++axis) {
                    frontWinding->p[frontWinding->numpoints][axis] =
                        splitPoint[axis];
                    backWinding->p[backWinding->numpoints][axis] =
                        splitPoint[axis];
                }
                frontWinding->numpoints++;
                backWinding->numpoints++;
            }
        }
    }

    if (splitCapacity < frontWinding->numpoints ||
        splitCapacity < backWinding->numpoints) {
        Com_Error(1, "\x15" "ClipWinding: points exceeded estimate");
    }
    if (CODUO_CM_WINDING_MAX_POINTS < frontWinding->numpoints ||
        CODUO_CM_WINDING_MAX_POINTS < backWinding->numpoints) {
        Com_Error(1, "\x15" "ClipWinding: MAX_POINTS_ON_WINDING");
    }
}

void ChopWindingInPlace(winding_t **winding, const vec3_t normal,
                        float dist, float epsilon)
{
    winding_t *source = *winding;
    int32_t sideCounts[3] = {0, 0, 0};
    int32_t sides[CODUO_CM_WINDING_TEMP_POINT_LIMIT];
    float dists[CODUO_CM_WINDING_TEMP_POINT_LIMIT];

    for (int32_t pointIndex = 0; pointIndex < source->numpoints;
         ++pointIndex) {
#if EMULATE_X87
        /* x87-faithful: 80-bit dot accumulate (x,y,z), round to the 32-bit
         * cm_windingChopDist slot (fstps), reload, subtract dist in 80-bit,
         * round again (flds;fsubs;fstps) — see the original instruction stream
         * and coduo_x87emu.h. */
        x87f cd = x87f_mul(x87f_load_f32(source->p[pointIndex][0]),
                           x87f_load_f32(normal[0]));
        cd = x87f_add(cd, x87f_mul(x87f_load_f32(source->p[pointIndex][1]),
                                   x87f_load_f32(normal[1])));
        cd = x87f_add(cd, x87f_mul(x87f_load_f32(source->p[pointIndex][2]),
                                   x87f_load_f32(normal[2])));
        cm_windingChopDist = x87f_store_f32(cd);
        cm_windingChopDist =
            x87f_store_f32(x87f_sub(x87f_load_f32(cm_windingChopDist),
                                    x87f_load_f32(dist)));
#else
        cm_windingChopDist =
            (source->p[pointIndex][0] * normal[0]) +
            (source->p[pointIndex][1] * normal[1]) +
            (source->p[pointIndex][2] * normal[2]);
        cm_windingChopDist -= dist;
#endif
        dists[pointIndex] = cm_windingChopDist;

        if (epsilon < cm_windingChopDist) {
            sides[pointIndex] = CODUO_CM_PATCH_PLANE_SIDE_FRONT;
        } else if (cm_windingChopDist < -epsilon) {
            sides[pointIndex] = CODUO_CM_PATCH_PLANE_SIDE_BACK;
        } else {
            sides[pointIndex] = CODUO_CM_PATCH_PLANE_SIDE_ON;
        }
        sideCounts[sides[pointIndex]]++;
    }

    sides[source->numpoints] = sides[0];
    dists[source->numpoints] = dists[0];

    if (sideCounts[CODUO_CM_PATCH_PLANE_SIDE_FRONT] == 0) {
        FreeWinding(source);
        *winding = NULL;
        return;
    }

    if (sideCounts[CODUO_CM_PATCH_PLANE_SIDE_BACK] == 0) {
        return;
    }

    int32_t splitCapacity =
        source->numpoints + CODUO_CM_WINDING_SPLIT_EXTRA_POINTS;
    winding_t *frontWinding = AllocWinding(splitCapacity);

    for (int32_t pointIndex = 0; pointIndex < source->numpoints;
         ++pointIndex) {
        const float *point = source->p[pointIndex];

        if (sides[pointIndex] == CODUO_CM_PATCH_PLANE_SIDE_ON) {
            for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT;
                 ++axis) {
                frontWinding->p[frontWinding->numpoints][axis] =
                    point[axis];
            }
            frontWinding->numpoints++;
        } else {
            if (sides[pointIndex] == CODUO_CM_PATCH_PLANE_SIDE_FRONT) {
                for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT;
                     ++axis) {
                    frontWinding->p[frontWinding->numpoints][axis] =
                        point[axis];
                }
                frontWinding->numpoints++;
            }

            if (sides[pointIndex + 1] != CODUO_CM_PATCH_PLANE_SIDE_ON &&
                sides[pointIndex + 1] != sides[pointIndex]) {
                const float *nextPoint =
                    source->p[(pointIndex + 1) % source->numpoints];
#if EMULATE_X87
                /* x87-faithful split ratio: dists[i] / (dists[i] -
                 * dists[i+1]), computed in 80-bit and rounded to the 32-bit
                 * cm_windingChopDist slot. */
                cm_windingChopDist = x87f_store_f32(
                    x87f_div(x87f_load_f32(dists[pointIndex]),
                             x87f_sub(x87f_load_f32(dists[pointIndex]),
                                      x87f_load_f32(dists[pointIndex + 1]))));
#else
                cm_windingChopDist =
                    dists[pointIndex] /
                    (dists[pointIndex] - dists[pointIndex + 1]);
#endif

                vec3_t splitPoint;
                for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT;
                     ++axis) {
                    if (normal[axis] == CODUO_CM_WINDING_AXIS_NORMAL) {
                        splitPoint[axis] = dist;
                    } else if (normal[axis] == -CODUO_CM_WINDING_AXIS_NORMAL) {
                        splitPoint[axis] = -dist;
                    } else {
#if EMULATE_X87
                        /* point + (nextPoint - point) * ratio, in 80-bit,
                         * rounded once to the 32-bit split vertex slot. */
                        splitPoint[axis] = x87f_store_f32(x87f_add(
                            x87f_load_f32(point[axis]),
                            x87f_mul(x87f_sub(x87f_load_f32(nextPoint[axis]),
                                              x87f_load_f32(point[axis])),
                                     x87f_load_f32(cm_windingChopDist))));
#else
                        splitPoint[axis] =
                            point[axis] +
                            (nextPoint[axis] - point[axis]) *
                                cm_windingChopDist;
#endif
                    }
                }

                for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT;
                     ++axis) {
                    frontWinding->p[frontWinding->numpoints][axis] =
                        splitPoint[axis];
                }
                frontWinding->numpoints++;
            }
        }
    }

    if (splitCapacity < frontWinding->numpoints) {
        Com_Error(1, "\x15" "ClipWinding: points exceeded estimate");
    }
    if (CODUO_CM_WINDING_MAX_POINTS < frontWinding->numpoints) {
        Com_Error(1, "\x15" "ClipWinding: MAX_POINTS_ON_WINDING");
    }

    FreeWinding(source);
    *winding = frontWinding;
}

winding_t *ChopWinding(winding_t *winding,
                                 const vec3_t normal, float dist)
{
    winding_t *front;
    winding_t *back;

    ClipWindingEpsilon(winding, normal, dist, CODUO_CM_WINDING_CLIP_EPSILON,
                       &front, &back);
    FreeWinding(winding);
    if (back != NULL) {
        FreeWinding(back);
    }

    return front;
}

void CheckWinding(const winding_t *winding)
{
    if (winding->numpoints < CODUO_CM_WINDING_MIN_POINT_COUNT) {
        Com_Error(1, "\x15" "CheckWinding: %i points", winding->numpoints);
    }

    float area = WindingArea(winding);
    if (area < CODUO_CM_WINDING_MIN_AREA) {
        Com_Error(1, "\x15" "CheckWinding: %f area", (double)area);
    }

    vec3_t normal;
    float dist;
    WindingPlane(winding, normal, &dist);

    for (int32_t pointIndex = 0; pointIndex < winding->numpoints;
         ++pointIndex) {
        const float *point = winding->p[pointIndex];
        for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT;
             ++axis) {
            if ((float)CODUO_CM_WINDING_WORLD_COORD_LIMIT < point[axis] ||
                point[axis] < (float)-CODUO_CM_WINDING_WORLD_COORD_LIMIT) {
                Com_Error(1, "\x15" "CheckFace: BUGUS_RANGE: %f",
                          (double)point[axis]);
            }
        }

        float planeDist = (point[0] * normal[0]) + (point[1] * normal[1]) +
                          (point[2] * normal[2]) - dist;
        if (planeDist < -CODUO_CM_WINDING_CLIP_EPSILON ||
            CODUO_CM_WINDING_CLIP_EPSILON < planeDist) {
            Com_Error(1, "\x15" "CheckWinding: point off plane");
        }

        int32_t nextIndex =
            pointIndex + 1 == winding->numpoints ? 0 : pointIndex + 1;
        const float *nextPoint = winding->p[nextIndex];
        vec3_t edge = {
            nextPoint[0] - point[0],
            nextPoint[1] - point[1],
            nextPoint[2] - point[2],
        };

        float edgeLength =
            (float)sqrt((double)((edge[0] * edge[0]) + (edge[1] * edge[1]) +
                                 (edge[2] * edge[2])));
        if (edgeLength < CODUO_CM_WINDING_CLIP_EPSILON) {
            Com_Error(1, "\x15" "CheckWinding: degenerate edge");
        }

        vec3_t edgeNormal;
        CrossProduct(normal, edge, edgeNormal);
        VectorNormalize2(edgeNormal, edgeNormal);
        /* The original rounds the dot to the edgeDist slot BEFORE adding the
         * epsilon (fstp 0x8052fd4; fld; fadd 0.1f; fstp 0x8052fe2) - two
         * roundings, so keep the += as a separate statement. */
        float edgeDist = (point[0] * edgeNormal[0]) +
                         (point[1] * edgeNormal[1]) +
                         (point[2] * edgeNormal[2]);
        edgeDist += CODUO_CM_WINDING_CLIP_EPSILON;

        for (int32_t checkIndex = 0; checkIndex < winding->numpoints;
             ++checkIndex) {
            if (checkIndex == pointIndex) {
                continue;
            }

            const float *checkPoint = winding->p[checkIndex];
            float checkDist = (checkPoint[0] * edgeNormal[0]) +
                              (checkPoint[1] * edgeNormal[1]) +
                              (checkPoint[2] * edgeNormal[2]);
            if (edgeDist < checkDist) {
                Com_Error(1, "\x15" "CheckWinding: non-convex");
            }
        }
    }
}

int32_t WindingOnPlaneSide(const winding_t *winding,
                           const vec3_t normal, float dist)
{
    qboolean foundFront = qfalse;
    qboolean foundBack = qfalse;

    for (int32_t pointIndex = 0; pointIndex < winding->numpoints;
         ++pointIndex) {
        float pointDist = (winding->p[pointIndex][0] * normal[0]) +
                          (winding->p[pointIndex][1] * normal[1]) +
                          (winding->p[pointIndex][2] * normal[2]) -
                          dist;

        if (pointDist < -CODUO_CM_WINDING_CLIP_EPSILON) {
            if (foundFront != qfalse) {
                return CODUO_CM_PATCH_PLANE_SIDE_CROSS;
            }
            foundBack = qtrue;
        } else if (CODUO_CM_WINDING_CLIP_EPSILON < pointDist) {
            if (foundBack != qfalse) {
                return CODUO_CM_PATCH_PLANE_SIDE_CROSS;
            }
            foundFront = qtrue;
        }
    }

    if (foundBack != qfalse) {
        return CODUO_CM_PATCH_PLANE_SIDE_BACK;
    }

    if (foundFront != qfalse) {
        return CODUO_CM_PATCH_PLANE_SIDE_FRONT;
    }

    return CODUO_CM_PATCH_PLANE_SIDE_ON;
}

void AddWindingToConvexHull(const winding_t *source,
                  winding_t **winding, const vec3_t normal)
{
    if (*winding == NULL) {
        *winding = CopyWinding(source);
        return;
    }

    int32_t mergedCount = (*winding)->numpoints;
    vec3_t mergedPoints[CODUO_CM_WINDING_MERGE_POINT_LIMIT];
    vec3_t newPoints[CODUO_CM_WINDING_MERGE_POINT_LIMIT];
    vec3_t edgeNormals[CODUO_CM_WINDING_MERGE_POINT_LIMIT];
    qboolean keepEdge[CODUO_CM_WINDING_MERGE_POINT_LIMIT];

    uint32_t mergedBytes =
        (uint32_t)mergedCount * (uint32_t)sizeof(mergedPoints[0]);
    memcpy(mergedPoints, (*winding)->p, (size_t)mergedBytes);

    for (int32_t sourceIndex = 0; sourceIndex < source->numpoints;
         ++sourceIndex) {
        const float *sourcePoint = source->p[sourceIndex];

        for (int32_t pointIndex = 0; pointIndex < mergedCount; ++pointIndex) {
            int32_t nextIndex = (pointIndex + 1) % mergedCount;
            vec3_t edge = {
                mergedPoints[nextIndex][0] - mergedPoints[pointIndex][0],
                mergedPoints[nextIndex][1] - mergedPoints[pointIndex][1],
                mergedPoints[nextIndex][2] - mergedPoints[pointIndex][2],
            };

            VectorNormalize2(edge, edge);
            CrossProduct(normal, edge, edgeNormals[pointIndex]);
        }

        qboolean pointOutside = qfalse;
        for (int32_t pointIndex = 0; pointIndex < mergedCount; ++pointIndex) {
            vec3_t delta = {
                sourcePoint[0] - mergedPoints[pointIndex][0],
                sourcePoint[1] - mergedPoints[pointIndex][1],
                sourcePoint[2] - mergedPoints[pointIndex][2],
            };
            float edgeDist = (delta[0] * edgeNormals[pointIndex][0]) +
                             (delta[1] * edgeNormals[pointIndex][1]) +
                             (delta[2] * edgeNormals[pointIndex][2]);

            if (CODUO_CM_WINDING_CLIP_EPSILON <= edgeDist) {
                pointOutside = qtrue;
            }
            keepEdge[pointIndex] =
                -CODUO_CM_WINDING_CLIP_EPSILON <= edgeDist ? qtrue : qfalse;
        }

        if (pointOutside == qfalse) {
            continue;
        }

        int32_t startEdge = 0;
        while (startEdge < mergedCount &&
               (keepEdge[startEdge % mergedCount] != qfalse ||
                keepEdge[(startEdge + 1) % mergedCount] == qfalse)) {
            startEdge++;
        }

        if (startEdge == mergedCount) {
            continue;
        }

        newPoints[0][0] = sourcePoint[0];
        newPoints[0][1] = sourcePoint[1];
        newPoints[0][2] = sourcePoint[2];
        int32_t newCount = 1;
        startEdge = (startEdge + 1) % mergedCount;

        for (int32_t pointOffset = 0; pointOffset < mergedCount;
             ++pointOffset) {
            int32_t edgeIndex = (startEdge + pointOffset) % mergedCount;
            int32_t nextIndex = (startEdge + pointOffset + 1) % mergedCount;

            if (keepEdge[edgeIndex] == qfalse ||
                keepEdge[nextIndex] == qfalse) {
                newPoints[newCount][0] = mergedPoints[nextIndex][0];
                newPoints[newCount][1] = mergedPoints[nextIndex][1];
                newPoints[newCount][2] = mergedPoints[nextIndex][2];
                newCount++;
            }
        }

        mergedCount = newCount;
        uint32_t newBytes =
            (uint32_t)newCount * (uint32_t)sizeof(newPoints[0]);
        memcpy(mergedPoints, newPoints, (size_t)newBytes);
    }

    FreeWinding(*winding);
    winding_t *merged = AllocWinding(mergedCount);
    merged->numpoints = mergedCount;
    *winding = merged;
    mergedBytes =
        (uint32_t)mergedCount * (uint32_t)sizeof(mergedPoints[0]);
    memcpy(merged->p, mergedPoints, (size_t)mergedBytes);
}

/* NOT_FROM_ORIGINAL_SOURCE: typed access to stored patch planes with inward flips. */
static void CM_LoadPatchFacetPlane(int32_t planeIndex, qboolean inward,
                                   vec4_t plane)
{
    plane[0] = cm_patchPlanes[planeIndex].normal[0];
    plane[1] = cm_patchPlanes[planeIndex].normal[1];
    plane[2] = cm_patchPlanes[planeIndex].normal[2];
    plane[3] = cm_patchPlanes[planeIndex].dist;

    if (inward == qfalse) {
        plane[0] = -plane[0];
        plane[1] = -plane[1];
        plane[2] = -plane[2];
        plane[3] = -plane[3];
    }
}

qboolean CM_ValidateFacet(const facet_t *facet)
{
    if (facet->surfacePlane == -1) {
        return qfalse;
    }

    vec4_t plane;
    CM_LoadPatchFacetPlane(facet->surfacePlane, qtrue, plane);
    winding_t *winding = BaseWindingForPlane(plane, plane[3]);

    for (int32_t borderIndex = 0;
         borderIndex < facet->numBorders && winding != NULL;
         ++borderIndex) {
        if (facet->borderPlanes[borderIndex] == -1) {
            FreeWinding(winding);
            return qfalse;
        }

        CM_LoadPatchFacetPlane(facet->borderPlanes[borderIndex],
                               facet->borderInward[borderIndex], plane);
        ChopWindingInPlace(&winding, plane, plane[3],
                           CODUO_CM_PATCH_BEVEL_PLANE_EPSILON);
    }

    if (winding == NULL) {
        return qfalse;
    }

    vec3_t mins;
    vec3_t maxs;
    WindingBounds(winding, mins, maxs);
    FreeWinding(winding);

    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        if ((maxs[axis] - mins[axis]) >
            (float)CODUO_CM_WINDING_WORLD_COORD_LIMIT) {
            return qfalse;
        }
        if ((float)CODUO_CM_WINDING_WORLD_COORD_LIMIT <= mins[axis]) {
            return qfalse;
        }
        if (maxs[axis] <= (float)-CODUO_CM_WINDING_WORLD_COORD_LIMIT) {
            return qfalse;
        }
    }

    return qtrue;
}

void CM_AddFacetBevels(facet_t *facet)
{
    vec4_t surfacePlane;
    CM_LoadPatchFacetPlane(facet->surfacePlane, qtrue, surfacePlane);
    winding_t *winding =
        BaseWindingForPlane(surfacePlane, surfacePlane[3]);

    for (int32_t borderIndex = 0;
         borderIndex < facet->numBorders && winding != NULL;
         ++borderIndex) {
        if (facet->borderPlanes[borderIndex] == facet->surfacePlane) {
            continue;
        }

        vec4_t plane;
        CM_LoadPatchFacetPlane(facet->borderPlanes[borderIndex],
                               facet->borderInward[borderIndex], plane);
        ChopWindingInPlace(&winding, plane, plane[3],
                           CODUO_CM_PATCH_BEVEL_PLANE_EPSILON);
    }

    if (winding == NULL) {
        return;
    }

    vec3_t mins;
    vec3_t maxs;
    WindingBounds(winding, mins, maxs);

    for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
        for (int32_t direction = CODUO_CM_PATCH_BEVEL_AXIS_MIN;
             direction <= CODUO_CM_PATCH_BEVEL_AXIS_MAX; direction += 2) {
            vec4_t bevel = {0.0f, 0.0f, 0.0f, 0.0f};
            bevel[axis] = (float)direction;
            bevel[3] = direction == 1 ? maxs[axis] : -mins[axis];

            qboolean flipped;
            if (CM_PlaneEqual(surfacePlane, bevel, &flipped) != qfalse) {
                continue;
            }

            int32_t borderIndex;
            for (borderIndex = 0; borderIndex < facet->numBorders;
                 ++borderIndex) {
                vec4_t border;
                CM_LoadPatchFacetPlane(facet->borderPlanes[borderIndex], qtrue,
                                       border);
                if (CM_PlaneEqual(border, bevel, &flipped) != qfalse) {
                    break;
                }
            }

            if (borderIndex != facet->numBorders) {
                continue;
            }

            /*
             * ORIGINAL_BINARY_BUG: the stock check is strictly `> 26`, only
             * prints, and then writes anyway. Index 26 already crosses the
             * fixed arrays. See
             * original-bugs/coduomp-collision-patch-facet-border-overflow.md.
             */
            if (facet->numBorders >
                CM_PATCH_FACET_BORDER_LIMIT) {
                Com_Printf("ERROR: too many bevels\n");
            }

            int32_t insert = facet->numBorders;
            facet->borderPlanes[insert] = CM_FindPlane(bevel, &flipped);
            facet->borderNoAdjust[insert] = qfalse;
            facet->borderInward[insert] = flipped;
            facet->numBorders++;
        }
    }

    for (int32_t pointIndex = 0; pointIndex < winding->numpoints;
         ++pointIndex) {
        int32_t nextIndex = (pointIndex + 1) % winding->numpoints;
#if EMULATE_X87
        /* edge = points[i] - points[next], each an 80-bit subtract rounded to
         * the float slot (DLL 0x0804e76b). */
        vec3_t edge;
        for (int32_t k = 0; k < CM_PATCH_VECTOR_AXIS_COUNT; ++k) {
            edge[k] = x87f_store_f32(
                x87f_sub(x87f_load_f32(winding->p[pointIndex][k]),
                         x87f_load_f32(winding->p[nextIndex][k])));
        }
#else
        vec3_t edge = {
            winding->p[pointIndex][0] - winding->p[nextIndex][0],
            winding->p[pointIndex][1] - winding->p[nextIndex][1],
            winding->p[pointIndex][2] - winding->p[nextIndex][2],
        };
#endif

        if (VectorNormalize(edge) < CODUO_CM_PATCH_BEVEL_EDGE_LENGTH_MIN) {
            continue;
        }

        CM_SnapVector(edge);

        int32_t axis;
        for (axis = 0;
             axis < CM_PATCH_VECTOR_AXIS_COUNT && edge[axis] != -1.0f &&
             edge[axis] != 1.0f;
             ++axis) {
        }
        if (axis < CM_PATCH_VECTOR_AXIS_COUNT) {
            continue;
        }

        for (axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
            for (int32_t direction = CODUO_CM_PATCH_BEVEL_AXIS_MIN;
                 direction <= CODUO_CM_PATCH_BEVEL_AXIS_MAX; direction += 2) {
                vec3_t axial = {0.0f, 0.0f, 0.0f};
                axial[axis] = (float)direction;

                vec4_t bevel;
                CrossProduct(edge, axial, bevel);
                if (VectorNormalize(bevel) <
                    CODUO_CM_PATCH_BEVEL_EDGE_LENGTH_MIN) {
                    continue;
                }

#if EMULATE_X87
                /* bevel[3] = winding->p[pointIndex] . bevel: three 80-bit
                 * products, two 80-bit adds, rounded to float once (DLL
                 * 0x0804e8f0). */
                bevel[3] = x87f_store_f32(x87f_add(
                    x87f_add(x87f_mul(x87f_load_f32(winding->p[pointIndex][0]),
                                      x87f_load_f32(bevel[0])),
                             x87f_mul(x87f_load_f32(winding->p[pointIndex][1]),
                                      x87f_load_f32(bevel[1]))),
                    x87f_mul(x87f_load_f32(winding->p[pointIndex][2]),
                             x87f_load_f32(bevel[2]))));
#else
                bevel[3] = (winding->p[pointIndex][0] * bevel[0]) +
                            (winding->p[pointIndex][1] * bevel[1]) +
                            (winding->p[pointIndex][2] * bevel[2]);
#endif

                qboolean hasBackPoint = qfalse;
                int32_t checkPoint;
                for (checkPoint = 0; checkPoint < winding->numpoints;
                     ++checkPoint) {
#if EMULATE_X87
                    /* dist = points[checkPoint] . bevel - bevel[3]: one 80-bit
                     * chain rounded to the float dist slot (DLL 0x0804e962),
                     * then compared against the double bevel-plane epsilon. */
                    float dist = x87f_store_f32(x87f_sub(
                        x87f_add(
                            x87f_add(
                                x87f_mul(x87f_load_f32(
                                             winding->p[checkPoint][0]),
                                         x87f_load_f32(bevel[0])),
                                x87f_mul(x87f_load_f32(
                                             winding->p[checkPoint][1]),
                                         x87f_load_f32(bevel[1]))),
                            x87f_mul(x87f_load_f32(
                                         winding->p[checkPoint][2]),
                                     x87f_load_f32(bevel[2]))),
                        x87f_load_f32(bevel[3])));
                    if (x87f_lt(x87f_load_f64(CODUO_CM_PATCH_BEVEL_PLANE_EPSILON),
                                x87f_load_f32(dist))) {
                        break;
                    }
                    if (x87f_lt(x87f_load_f32(dist),
                                x87f_load_f64(-CODUO_CM_PATCH_BEVEL_PLANE_EPSILON))) {
                        hasBackPoint = qtrue;
                    }
#else
                    float dist =
                        (winding->p[checkPoint][0] * bevel[0]) +
                        (winding->p[checkPoint][1] * bevel[1]) +
                        (winding->p[checkPoint][2] * bevel[2]) -
                        bevel[3];
                    if (dist > CODUO_CM_PATCH_BEVEL_PLANE_EPSILON) {
                        break;
                    }
                    if (dist < -CODUO_CM_PATCH_BEVEL_PLANE_EPSILON) {
                        hasBackPoint = qtrue;
                    }
#endif
                }

                if (checkPoint != winding->numpoints ||
                    hasBackPoint == qfalse) {
                    continue;
                }

                int32_t borderIndex;
                qboolean flipped;
                for (borderIndex = 0; borderIndex < facet->numBorders;
                     ++borderIndex) {
                    vec4_t border;
                    CM_LoadPatchFacetPlane(facet->borderPlanes[borderIndex],
                                           qtrue, border);
                    if (CM_PlaneEqual(border, bevel, &flipped) != qfalse) {
                        break;
                    }
                }

                if (borderIndex != facet->numBorders) {
                    continue;
                }
                /*
                 * ORIGINAL_BINARY_BUG: this is the second identical late,
                 * non-blocking stock check; the insertion still follows.
                 */
                if (facet->numBorders >
                    CM_PATCH_FACET_BORDER_LIMIT) {
                    Com_Printf("ERROR: too many bevels\n");
                }

                int32_t insert = facet->numBorders;
                facet->borderPlanes[insert] = CM_FindPlane(bevel, &flipped);
                for (int32_t prior = 0; prior < facet->numBorders;
                     ++prior) {
                    if (facet->borderPlanes[insert] ==
                        facet->borderPlanes[prior]) {
                        Com_Printf("WARNING: bevel plane already used\n");
                    }
                }

                facet->borderNoAdjust[insert] = qfalse;
                facet->borderInward[insert] = flipped;

                winding_t *test = CopyWinding(winding);
                vec4_t plane;
                CM_LoadPatchFacetPlane(facet->borderPlanes[insert],
                                       facet->borderInward[insert], plane);
                ChopWindingInPlace(&test, plane, plane[3],
                                   CODUO_CM_PATCH_BEVEL_PLANE_EPSILON);
                if (test == NULL) {
                    Com_DPrintf(
                        "WARNING: CM_AddFacetBevels... invalid bevel\n");
                } else {
                    FreeWinding(test);
                    facet->numBorders++;
                }
            }
        }
    }

    FreeWinding(winding);

    /*
     * ORIGINAL_BINARY_BUG: stock appends the surface plane without any
     * capacity check, even after the non-blocking bevel warnings above.
     */
    int32_t insert = facet->numBorders;
    facet->borderPlanes[insert] = facet->surfacePlane;
    facet->borderNoAdjust[insert] = qfalse;
    facet->borderInward[insert] = qtrue;
    facet->numBorders++;
}

void CM_PatchCollideFromGrid(const cGrid_t *grid,
                             patchCollide_t *patchCollide)
{
    cm_patchPlaneCount = 0;
    int32_t numFacets = 0;

    for (int32_t x = 0; x < grid->width - 1; ++x) {
        for (int32_t y = 0; y < grid->height - 1; ++y) {
            cm_patchPlaneGrid[x][y][0] =
                CM_FindPlane2(grid->points[x][y], grid->points[x + 1][y],
                              grid->points[x + 1][y + 1]);
            cm_patchPlaneGrid[x][y][1] =
                CM_FindPlane2(grid->points[x + 1][y + 1],
                              grid->points[x][y + 1], grid->points[x][y]);
        }
    }

    for (int32_t x = 0; x < grid->width - 1; ++x) {
        for (int32_t y = 0; y < grid->height - 1; ++y) {
            int32_t bottomPlane = -1;
            if (y == 0) {
                if (grid->wrapHeight != qfalse) {
                    bottomPlane = cm_patchPlaneGrid[x][grid->height - 2][1];
                }
            } else {
                bottomPlane = cm_patchPlaneGrid[x][y - 1][1];
            }
            qboolean bottomReused =
                bottomPlane == cm_patchPlaneGrid[x][y][0] ? qtrue : qfalse;
            if (bottomPlane == -1 || bottomReused != qfalse) {
                bottomPlane =
                    CM_EdgePlaneNum(grid, cm_patchPlaneGrid, x, y,
                                    CODUO_CM_PATCH_EDGE_BOTTOM);
            }

            int32_t topPlane = -1;
            if (y < grid->height - 2) {
                topPlane = cm_patchPlaneGrid[x][y + 1][0];
            } else if (grid->wrapHeight != qfalse) {
                topPlane = cm_patchPlaneGrid[x][0][0];
            }
            qboolean topReused =
                topPlane == cm_patchPlaneGrid[x][y][1] ? qtrue : qfalse;
            if (topPlane == -1 || topReused != qfalse) {
                topPlane = CM_EdgePlaneNum(grid, cm_patchPlaneGrid, x, y,
                                           CODUO_CM_PATCH_EDGE_TOP);
            }

            int32_t leftPlane = -1;
            if (x == 0) {
                if (grid->wrapWidth != qfalse) {
                    leftPlane = cm_patchPlaneGrid[grid->width - 2][y][0];
                }
            } else {
                leftPlane = cm_patchPlaneGrid[x - 1][y][0];
            }
            qboolean leftReused =
                leftPlane == cm_patchPlaneGrid[x][y][1] ? qtrue : qfalse;
            if (leftPlane == -1 || leftReused != qfalse) {
                leftPlane = CM_EdgePlaneNum(grid, cm_patchPlaneGrid, x, y,
                                            CODUO_CM_PATCH_EDGE_LEFT);
            }

            int32_t rightPlane = -1;
            if (x < grid->width - 2) {
                rightPlane = cm_patchPlaneGrid[x + 1][y][1];
            } else if (grid->wrapWidth != qfalse) {
                rightPlane = cm_patchPlaneGrid[0][y][1];
            }
            qboolean rightReused =
                rightPlane == cm_patchPlaneGrid[x][y][0] ? qtrue : qfalse;
            if (rightPlane == -1 || rightReused != qfalse) {
                rightPlane = CM_EdgePlaneNum(grid, cm_patchPlaneGrid, x, y,
                                             CODUO_CM_PATCH_EDGE_RIGHT);
            }

            if (numFacets == CODUO_CM_PATCH_MAX_FACETS) {
                Com_Error(1, "\x15" "MAX_FACETS");
            }

            facet_t *facet = &cm_patchFacets[numFacets];
            Com_Memset(facet, 0, sizeof(*facet));

            if (cm_patchPlaneGrid[x][y][0] == cm_patchPlaneGrid[x][y][1]) {
                if (cm_patchPlaneGrid[x][y][0] != -1) {
                    facet->surfacePlane = cm_patchPlaneGrid[x][y][0];
                    facet->numBorders = 4;
                    facet->borderPlanes[0] = bottomPlane;
                    facet->borderNoAdjust[0] = bottomReused;
                    facet->borderPlanes[1] = rightPlane;
                    facet->borderNoAdjust[1] = rightReused;
                    facet->borderPlanes[2] = topPlane;
                    facet->borderNoAdjust[2] = topReused;
                    facet->borderPlanes[3] = leftPlane;
                    facet->borderNoAdjust[3] = leftReused;
                    CM_SetBorderInward(facet, grid, cm_patchPlaneGrid, x,
                                       y, CODUO_CM_PATCH_BORDER_INWARD_QUAD);
                    if (CM_ValidateFacet(facet) != qfalse) {
                        CM_AddFacetBevels(facet);
                        numFacets++;
                    }
                }
                continue;
            }

            facet->surfacePlane = cm_patchPlaneGrid[x][y][0];
            facet->numBorders = 3;
            facet->borderPlanes[0] = bottomPlane;
            facet->borderNoAdjust[0] = bottomReused;
            facet->borderPlanes[1] = rightPlane;
            facet->borderNoAdjust[1] = rightReused;
            facet->borderPlanes[2] = cm_patchPlaneGrid[x][y][1];
            if (facet->borderPlanes[2] == -1) {
                facet->borderPlanes[2] =
                    topPlane != -1
                        ? topPlane
                        : CM_EdgePlaneNum(
                              grid, cm_patchPlaneGrid, x, y,
                              CODUO_CM_PATCH_EDGE_DIAGONAL_DESCENDING);
            }
            CM_SetBorderInward(facet, grid, cm_patchPlaneGrid, x, y,
                                CODUO_CM_PATCH_BORDER_INWARD_LOWER_TRIANGLE);
            if (CM_ValidateFacet(facet) != qfalse) {
                CM_AddFacetBevels(facet);
                numFacets++;
            }

            if (numFacets == CODUO_CM_PATCH_MAX_FACETS) {
                Com_Error(1, "\x15" "MAX_FACETS");
            }

            facet = &cm_patchFacets[numFacets];
            Com_Memset(facet, 0, sizeof(*facet));
            facet->surfacePlane = cm_patchPlaneGrid[x][y][1];
            facet->numBorders = 3;
            facet->borderPlanes[0] = topPlane;
            facet->borderNoAdjust[0] = topReused;
            facet->borderPlanes[1] = leftPlane;
            facet->borderNoAdjust[1] = leftReused;
            facet->borderPlanes[2] = cm_patchPlaneGrid[x][y][0];
            if (facet->borderPlanes[2] == -1) {
                facet->borderPlanes[2] =
                    bottomPlane != -1
                        ? bottomPlane
                        : CM_EdgePlaneNum(
                              grid, cm_patchPlaneGrid, x, y,
                              CODUO_CM_PATCH_EDGE_DIAGONAL_ASCENDING);
            }
            CM_SetBorderInward(facet, grid, cm_patchPlaneGrid, x, y,
                                CODUO_CM_PATCH_BORDER_INWARD_UPPER_TRIANGLE);
            if (CM_ValidateFacet(facet) != qfalse) {
                CM_AddFacetBevels(facet);
                numFacets++;
            }
        }
    }

    patchCollide->numPlanes = cm_patchPlaneCount;
    patchCollide->numFacets = numFacets;
    patchCollide->facets =
        Hunk_Alloc((size_t)numFacets * sizeof(patchCollide->facets[0]));
    Com_Memcpy(patchCollide->facets, cm_patchFacets,
                 (size_t)numFacets * sizeof(patchCollide->facets[0]));
    patchCollide->planes =
        Hunk_Alloc((size_t)cm_patchPlaneCount * sizeof(patchCollide->planes[0]));
    Com_Memcpy(patchCollide->planes, cm_patchPlanes,
                 (size_t)cm_patchPlaneCount * sizeof(patchCollide->planes[0]));
}

patchCollide_t *
CM_GeneratePatchCollide(uint32_t width, uint32_t height, int32_t maxError,
                        const vec3_t *points, vec3_t *bounds)
{
    if ((int32_t)width < CODUO_CM_WINDING_MIN_POINT_COUNT ||
        (int32_t)height < CODUO_CM_WINDING_MIN_POINT_COUNT ||
        points == NULL) {
        Com_Error(1, "\x15" "CM_GeneratePatchFacets: bad parameters: (%i, %i, %p)",
                  width, height, (const void *)points);
    }
    if ((width & 1U) == 0 || (height & 1U) == 0) {
        Com_Error(1,
                  "\x15" "CM_GeneratePatchFacets: even sizes are invalid for quadratic meshes");
    }
    if (CODUO_CM_PATCH_POINT_GRID_SIZE < (int32_t)width ||
        CODUO_CM_PATCH_POINT_GRID_SIZE < (int32_t)height) {
        Com_Error(1, "\x15" "CM_GeneratePatchFacets: source is > MAX_GRID_SIZE");
    }

    cGrid_t *grid = &cm_patchWorkGrid;
    grid->width = (int32_t)width;
    grid->height = (int32_t)height;
    grid->wrapWidth = qfalse;
    grid->wrapHeight = qfalse;

    for (int32_t x = 0; x < (int32_t)width; ++x) {
        for (int32_t y = 0; y < (int32_t)height; ++y) {
            grid->points[x][y][0] = points[y * width + x][0];
            grid->points[x][y][1] = points[y * width + x][1];
            grid->points[x][y][2] = points[y * width + x][2];
        }
    }

    CM_SetGridWrapWidth(grid);
    CM_SubdivideGridColumns(grid, maxError);
    CM_RemoveDegenerateColumns(grid);
    CM_TransposeGrid(grid);
    CM_SetGridWrapWidth(grid);
    CM_SubdivideGridColumns(grid, maxError);
    CM_RemoveDegenerateColumns(grid);

    patchCollide_t *patchCollide =
        Hunk_Alloc((size_t)sizeof(*patchCollide));

    ClearBounds(bounds[0], bounds[1]);
    for (int32_t x = 0; x < grid->width; ++x) {
        for (int32_t y = 0; y < grid->height; ++y) {
            AddPointToBounds(grid->points[x][y], bounds[0], bounds[1]);
        }
    }

    CM_PatchCollideFromGrid(grid, patchCollide);

    bounds[0][0] -= 1.0f;
    bounds[0][1] -= 1.0f;
    bounds[0][2] -= 1.0f;
    bounds[1][0] += 1.0f;
    bounds[1][1] += 1.0f;
    bounds[1][2] += 1.0f;

    return patchCollide;
}

#ifdef CODUO_COLLISION_DIGEST
/*
 * Curved-patch (terrainCollide) parity digest — NOT part of the reconstruction.
 *
 * The PATCHPARITY dump only covers patch->patchCollide (the triangle-soup
 * builder). Curved patches instead populate patch->terrainCollide via
 * CM_GeneratePatchCollide, whose float pipeline (CM_Subdivide, CM_PlaneFromPoints,
 * CM_FindPlane2, CM_AddFacetBevels, ChopWindingInPlace, ...) is otherwise
 * unmeasured. This hashes what that pipeline actually produces: every generated
 * plane's float bits, plus each facet's plane indices/flags (which catch
 * plane-selection and border-classification divergence).
 *
 * The build struct is private to this file, so the digest calls in here.
 */
void coduo_engine_collision_digest_bytes_external(const void *data, size_t length, uint64_t *accum);

/* NOT_FROM_ORIGINAL_SOURCE: curved-patch parity digest. */
void coduo_engine_digest_terrain_collide(const void *terrainCollide, uint64_t *accum,
                             int32_t *planeCountOut, int32_t *facetCountOut)
{
    const patchCollide_t *build = terrainCollide;

    *planeCountOut = build->numPlanes;
    *facetCountOut = build->numFacets;

    for (int32_t planeIndex = 0; planeIndex < build->numPlanes; ++planeIndex) {
        const patchPlane_t *plane = &build->planes[planeIndex];
        coduo_engine_collision_digest_bytes_external(plane->normal, sizeof(plane->normal), accum);
        coduo_engine_collision_digest_bytes_external(&plane->dist, sizeof(plane->dist), accum);
        coduo_engine_collision_digest_bytes_external(&plane->signbits, sizeof(plane->signbits),
                               accum);
    }

    for (int32_t facetIndex = 0; facetIndex < build->numFacets; ++facetIndex) {
        const facet_t *facet = &build->facets[facetIndex];
        coduo_engine_collision_digest_bytes_external(&facet->surfacePlane, sizeof(facet->surfacePlane),
                               accum);
        coduo_engine_collision_digest_bytes_external(&facet->numBorders,
                               sizeof(facet->numBorders), accum);
        for (int32_t border = 0; border < facet->numBorders; ++border) {
            coduo_engine_collision_digest_bytes_external(&facet->borderPlanes[border],
                                   sizeof(facet->borderPlanes[border]), accum);
            coduo_engine_collision_digest_bytes_external(&facet->borderInward[border],
                                   sizeof(facet->borderInward[border]), accum);
            coduo_engine_collision_digest_bytes_external(&facet->borderNoAdjust[border],
                                   sizeof(facet->borderNoAdjust[border]), accum);
        }
    }
}
#endif

void *CM_Hunk_AllocXModelMesh(size_t size)
{
    return Hunk_Alloc(size);
}

void *CM_Hunk_AllocXModel(size_t size)
{
    return Hunk_Alloc(size);
}

void CM_CreateStaticModel(worldSectorAreaLink_t *areaLink,
                          const char *modelName,
                          const vec3_t origin,
                          const vec3_t angles,
                          const vec3_t scale)
{
    if (modelName == NULL || modelName[0] == '\0') {
        Com_Error(1, "\x15" "Invalid static model name @ %f %f %f\n",
                  (double)origin[0], (double)origin[1],
                  (double)origin[2]);
    }
    if (scale[0] == 0.0f) {
        Com_Error(1, "\x15" "Static model [%s] has x scale of 0.0\n",
                  modelName);
    }
    if (scale[1] == 0.0f) {
        Com_Error(1, "\x15" "Static model [%s] has y scale of 0.0\n",
                  modelName);
    }
    if (scale[2] == 0.0f) {
        Com_Error(1, "\x15" "Static model [%s] has z scale of 0.0\n",
                  modelName);
    }

    areaLink->model = XModelPrecache(modelName, qtrue,
                                     CM_Hunk_AllocXModelMesh,
                                     CM_Hunk_AllocXModel);
    areaLink->origin[0] = origin[0];
    areaLink->origin[1] = origin[1];
    areaLink->origin[2] = origin[2];

    coduo_compact_affine_matrix43_t axis;
    AnglesToAxis(angles, &axis);

    vec3_t scaledAxis[3];
    for (int32_t row = 0; row < CM_PATCH_VECTOR_AXIS_COUNT; ++row) {
        for (int32_t column = 0; column < CM_PATCH_VECTOR_AXIS_COUNT;
             ++column) {
            scaledAxis[row][column] = axis.axis[row][column] * scale[row];
        }
    }

    MatrixInverse3x3(scaledAxis, areaLink->inverseAxis);

    if (XModelGetStaticBounds(areaLink->model, scaledAxis, areaLink->linkMins,
                              areaLink->linkMaxs) != qfalse) {
        for (int32_t axisIndex = 0; axisIndex < CM_PATCH_VECTOR_AXIS_COUNT;
             ++axisIndex) {
            areaLink->linkMins[axisIndex] += origin[axisIndex];
            areaLink->linkMaxs[axisIndex] += origin[axisIndex];
        }

        if (XModelGetContents(areaLink->model) != 0) {
            CM_LinkStaticModel(areaLink);
        }
    }
}

void CM_LoadStaticModels(void)
{
    char *parseData = cm_entityString;
    cm_staticModelCount = 0;
    cm_staticModels = NULL;

    for (;;) {
        const char *token = Com_Parse(&parseData);
        if (token[0] == '\0' || token[0] != '{') {
            break;
        }

        qboolean isMiscModel = qfalse;
        while ((token = Com_Parse(&parseData))[0] != '\0' &&
               token[0] != '}') {
            char key[CODUO_CM_MISC_MODEL_KEY_BUFFER_SIZE];
            char value[CODUO_CM_MISC_MODEL_VALUE_BUFFER_SIZE];

            /*
             * ORIGINAL_BINARY_BUG: stock copies parser tokens into these
             * fixed stack arrays without capacity checks. See
             * original-bugs/coduomp-static-model-entity-token-stack-overflow.md.
             */
            strcpy(key, token);
            token = Com_Parse(&parseData);
            if (token[0] == '\0') {
                break;
            }
            strcpy(value, token);

            if (strcasecmp(key, "classname") == 0 &&
                strcasecmp(value, "misc_model") == 0) {
                isMiscModel = qtrue;
            }
        }

        if (isMiscModel != qfalse) {
            cm_staticModelCount++;
        }
    }

    if (cm_staticModelCount == 0) {
        return;
    }

    cm_staticModels =
        Hunk_Alloc((size_t)cm_staticModelCount * sizeof(cm_staticModels[0]));
    parseData = cm_entityString;

    int32_t staticModelIndex = 0;
    for (;;) {
        const char *open = Com_Parse(&parseData);
        if (open[0] == '\0' || open[0] != '{') {
            break;
        }

        vec3_t origin = {0.0f, 0.0f, 0.0f};
        vec3_t angles = {0.0f, 0.0f, 0.0f};
        vec3_t scale = {1.0f, 1.0f, 1.0f};
        char modelName[CODUO_CM_MISC_MODEL_NAME_BUFFER_SIZE];
        qboolean isMiscModel = qfalse;

        modelName[0] = '\0';

        const char *token;
        while ((token = Com_Parse(&parseData))[0] != '\0' &&
               token[0] != '}') {
            char key[CODUO_CM_MISC_MODEL_KEY_BUFFER_SIZE];
            char value[CODUO_CM_MISC_MODEL_VALUE_BUFFER_SIZE];

            /*
             * ORIGINAL_BINARY_BUG: the second pass repeats the unbounded
             * fixed-array token copies documented in
             * original-bugs/coduomp-static-model-entity-token-stack-overflow.md.
             */
            strcpy(key, token);
            token = Com_Parse(&parseData);
            if (token[0] == '\0') {
                break;
            }
            strcpy(value, token);

            if (strcasecmp(key, "classname") == 0) {
                if (strcasecmp(value, "misc_model") == 0) {
                    isMiscModel = qtrue;
                }
            } else if (strcasecmp(key, "model") == 0) {
                /*
                 * ORIGINAL_BINARY_BUG: stock also copies the model token into
                 * a fixed 64-byte stack array without a bound. See
                 * original-bugs/coduomp-static-model-entity-token-stack-overflow.md.
                 */
                strcpy(modelName, value);
            } else if (strcasecmp(key, "origin") == 0) {
                sscanf(value, "%f %f %f", &origin[0], &origin[1],
                       &origin[2]);
            } else if (strcasecmp(key, "angles") == 0) {
                sscanf(value, "%f %f %f", &angles[0], &angles[1],
                       &angles[2]);
            } else if (strcasecmp(key, "modelscale_vec") == 0) {
                sscanf(value, "%f %f %f", &scale[0], &scale[1], &scale[2]);
            } else if (strcasecmp(key, "modelscale") == 0) {
                scale[0] = (float)atof(value);
                scale[1] = scale[0];
                scale[2] = scale[0];
            }
        }

        if (isMiscModel != qfalse) {
            /*
             * ORIGINAL_BINARY_BUG: stock skips seven bytes without checking
             * the value length or an xmodel/ prefix. See
             * original-bugs/coduomp-static-model-entity-token-stack-overflow.md.
             */
            CM_CreateStaticModel(&cm_staticModels[staticModelIndex],
                                 modelName + 7, origin, angles, scale);
            staticModelIndex++;
        }
    }
}

void CM_TraceStaticModel(worldSectorAreaLink_t *areaLink,
                         trace_t *trace,
                         const vec3_t start,
                         const vec3_t end,
                         int32_t contentsMask)
{
    /*
     * ORIGINAL_BINARY_BUG: original-machine-code parity depends on this partial
     * initialization. CM_TraceStaticModel at 0x8053da6..0x8053dab seeds only
     * the local trace fraction immediately on entry, then copies the full
     * 0x30-byte local trace record back at 0x8053f89..0x8053f99. Under the
     * current trace layout this leaves material, partName, and partGroup
     * metadata dependent on stack contents for static-model hits. See
     * original-bugs/coduomp-static-model-trace-uninitialized-metadata.md.
     * Do not initialize those fields in the recovered stock path.
     */
    trace_t localTrace;
    localTrace.fraction = trace->fraction;

    int32_t partNameCount =
        XModelNumBones(areaLink->model);
    DObjSkelMat basePose[partNameCount];
    XModelGetBasePose(areaLink->model, basePose);

    vec3_t offsetStart = {
        start[0] - areaLink->origin[0],
        start[1] - areaLink->origin[1],
        start[2] - areaLink->origin[2],
    };
    vec3_t localStart;
    MatrixTransformVector(offsetStart, areaLink->inverseAxis, localStart);

    vec3_t offsetEnd = {
        end[0] - areaLink->origin[0],
        end[1] - areaLink->origin[1],
        end[2] - areaLink->origin[2],
    };
    vec3_t localEnd;
    MatrixTransformVector(offsetEnd, areaLink->inverseAxis, localEnd);

    if (XModelTraceLine(areaLink->model, &localTrace, basePose, localStart,
                        localEnd, contentsMask) >= 0) {
        localTrace.entityNum = CODUO_CM_ENTITYNUM_WORLD;

        /* The original rounds each end-start delta to a float slot
         * (0x8053ecd/0x8053ee3/0x8053ef9) before the lerp multiply; keep the
         * explicit delta temporaries so the lerp reloads the rounded floats. */
        vec3_t delta = {
            end[0] - start[0],
            end[1] - start[1],
            end[2] - start[2],
        };
        localTrace.endpos[0] = start[0] + delta[0] * localTrace.fraction;
        localTrace.endpos[1] = start[1] + delta[1] * localTrace.fraction;
        localTrace.endpos[2] = start[2] + delta[2] * localTrace.fraction;

        vec3_t worldNormal;
        MatrixTransposeTransformVector(localTrace.normal, areaLink->inverseAxis,
                                       worldNormal);
        VectorNormalize(worldNormal);
        localTrace.normal[0] = worldNormal[0];
        localTrace.normal[1] = worldNormal[1];
        localTrace.normal[2] = worldNormal[2];

        *trace = localTrace;
    }
}

int32_t FindEdge(int32_t edgeCount,
                 const coduo_cm_triangle_edge_t *edges,
                 int16_t first,
                 int16_t second)
{
    for (int32_t edgeIndex = 0; edgeIndex < edgeCount; ++edgeIndex) {
        if ((edges[edgeIndex].first == second &&
             edges[edgeIndex].second == first) ||
            (edges[edgeIndex].first == first &&
             edges[edgeIndex].second == second)) {
            return edgeIndex;
        }
    }

    return -1;
}

coduo_patch_collide_t *
CM_GenerateTerrainCollide(int32_t indexCount, const int16_t *indices,
                          uint32_t vertexCount, const vec3_t *vertices,
                          vec3_t *bounds)
{
    if (indexCount % CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT != 0) {
        Com_Error(1,
                  "\x15" "CM_GenerateTerrainCollide: numIndexes % 3 != 0, corrupt bsp?");
    }
    if (CODUO_CM_MODEL_VERTEX_REMAP_LIMIT - 1U < vertexCount) {
        Com_Error(1, "\x15" "CM_GenerateTerrainCollide: too many vertices");
    }

    uint32_t edgeCount = 0;
    /*
     * ORIGINAL_BINARY_BUG: signed authored indices are used directly as
     * vertex and remap-table subscripts without a vertex-count check. See
     * original-bugs/coduomp-collision-bsp-cross-reference-bounds.md.
     */
    for (int32_t triangleIndex = 0; triangleIndex < indexCount;
         triangleIndex += CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT) {
        for (int32_t corner = 0; corner < CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT;
             ++corner) {
            int16_t first = indices[triangleIndex + corner];
            int16_t second =
                indices[triangleIndex +
                        ((corner + 1) % CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT)];
            int16_t opposite =
                indices[triangleIndex +
                        ((corner + 2) % CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT)];

            if (first == second) {
                Com_Error(1,
                          "\x15" "CM_GenerateTerrainCollide: degenerate triangle, corrupt bsp?");
            }

            int32_t edgeIndex =
                FindEdge((int32_t)edgeCount, cm_modelEdges, first, second);
            if (edgeIndex < 0) {
                if (CODUO_CM_MODEL_EDGE_LIMIT - 1U < edgeCount) {
                    Com_Error(1,
                              "\x15" "CM_GenerateTerrainCollide: too many edges");
                }
                cm_modelEdges[edgeCount].first = first;
                cm_modelEdges[edgeCount].second = second;
                cm_modelEdgeOppositeVertices[edgeCount] = opposite;
                edgeCount++;
                continue;
            }

            vec3_t firstDelta = {
                vertices[first][0] - vertices[opposite][0],
                vertices[first][1] - vertices[opposite][1],
                vertices[first][2] - vertices[opposite][2],
            };
            vec3_t secondDelta = {
                vertices[second][0] - vertices[opposite][0],
                vertices[second][1] - vertices[opposite][1],
                vertices[second][2] - vertices[opposite][2],
            };
            int16_t previousOpposite = cm_modelEdgeOppositeVertices[edgeIndex];
            vec3_t previousDelta = {
                vertices[previousOpposite][0] - vertices[opposite][0],
                vertices[previousOpposite][1] - vertices[opposite][1],
                vertices[previousOpposite][2] - vertices[opposite][2],
            };

            vec3_t cross;
            CrossProduct(secondDelta, firstDelta, cross);
            float concavity = (cross[0] * previousDelta[0]) +
                              (cross[1] * previousDelta[1]) +
                              (cross[2] * previousDelta[2]);
            if (concavity <= 0.0f) {
                continue;
            }

            float firstDistance =
                (float)DistanceSquared(vertices[first],
                                    vertices[previousOpposite]);
            float secondDistance =
                (float)DistanceSquared(vertices[second],
                                    vertices[previousOpposite]);
            float maxDistance =
                secondDistance < firstDistance ? firstDistance : secondDistance;
            /* The original never stores the cross-length dot: it is compared
             * in 80-bit straight against maxDistance*epsilon (0x805452a..
             * 0x8054544, no fstp), so it must stay inline in the condition. */
            if (((cross[0] * cross[0]) + (cross[1] * cross[1]) +
                 (cross[2] * cross[2])) <
                maxDistance * CODUO_CM_MODEL_COPLANAR_EDGE_EPSILON) {
                edgeCount--;
                cm_modelEdges[edgeIndex] = cm_modelEdges[edgeCount];
                cm_modelEdgeOppositeVertices[edgeIndex] =
                    cm_modelEdgeOppositeVertices[edgeCount];
            }
        }
    }

    memset(cm_modelVertexRemap, -1, sizeof(cm_modelVertexRemap));
    int32_t usedVertexCount = 0;
    for (int32_t edgeIndex = 0; edgeIndex < (int32_t)edgeCount; ++edgeIndex) {
        int16_t first = cm_modelEdges[edgeIndex].first;
        int16_t second = cm_modelEdges[edgeIndex].second;

        if (cm_modelVertexRemap[first] < 0) {
            cm_modelVertexRemap[first] = (int16_t)usedVertexCount;
            usedVertexCount++;
        }
        if (cm_modelVertexRemap[second] < 0) {
            cm_modelVertexRemap[second] = (int16_t)usedVertexCount;
            usedVertexCount++;
        }
    }

    int32_t triangleCount = indexCount / CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT;
    uint32_t soupBytes =
        (uint32_t)sizeof(coduo_cm_collision_triangle_soup_t) +
        (uint32_t)triangleCount * (uint32_t)sizeof(coduo_cm_collision_triangle_t);
    coduo_cm_collision_triangle_soup_t *soup =
        Hunk_Alloc(soupBytes);
    coduo_cm_collision_vertex_t *collisionVertices =
        Hunk_Alloc((size_t)usedVertexCount * sizeof(collisionVertices[0]));
    coduo_cm_collision_edge_t *collisionEdges =
        Hunk_Alloc((size_t)edgeCount * sizeof(collisionEdges[0]));

    qboolean hasPositiveZPlane = qfalse;
    qboolean hasNegativeZPlane = qfalse;
    soup->triangleCount = (uint16_t)triangleCount;
    const int32_t publishedTriangleCount = (int32_t)soup->triangleCount;

    for (int32_t triangleIndex = 0; triangleIndex < publishedTriangleCount;
         ++triangleIndex) {
        int16_t first =
            indices[triangleIndex * CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT];
        int16_t second =
            indices[triangleIndex * CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT + 1];
        int16_t third =
            indices[triangleIndex * CODUO_CM_MODEL_TRIANGLE_VERTEX_COUNT + 2];
        coduo_cm_collision_triangle_t *triangle =
            &soup->triangles[triangleIndex];

        PlaneFromPoints(triangle->plane, vertices[first], vertices[second],
                     vertices[third]);
        if (triangle->plane[2] < CODUO_CM_MODEL_TRIANGLE_PLANE_NEGATIVE_EPSILON) {
            hasNegativeZPlane = qtrue;
        } else if (CODUO_CM_MODEL_TRIANGLE_PLANE_POSITIVE_EPSILON <
                   triangle->plane[2]) {
            hasPositiveZPlane = qtrue;
        }

#if EMULATE_X87
        /* x87-faithful transcription of the svec/tvec basis (DLL 0x0805403e):
         * every product/sum/difference is an 80-bit chain rounded to its float
         * slot. VectorNormalize is the emulated variant; edge subtracts round
         * to float; dot is an 80-bit 3-term dot; inverseDot = 1/(1-dot*dot);
         * edge01/edge02 are scaled by inverseDot in place BEFORE svec/tvec;
         * each svec/tvec component is (edgeA - dot*edgeB) * (1/edgeLength);
         * svec[3]/tvec[3] are dots with vertices[first]. */
        vec3_t edge01;
        vec3_t edge02;
        for (int32_t k = 0; k < CM_PATCH_VECTOR_AXIS_COUNT; ++k) {
            edge01[k] = x87f_store_f32(x87f_sub(
                x87f_load_f32(vertices[second][k]),
                x87f_load_f32(vertices[first][k])));
        }
        float edge01Length = (float)VectorNormalize(edge01);
        for (int32_t k = 0; k < CM_PATCH_VECTOR_AXIS_COUNT; ++k) {
            edge02[k] = x87f_store_f32(x87f_sub(
                x87f_load_f32(vertices[third][k]),
                x87f_load_f32(vertices[first][k])));
        }
        float edge02Length = (float)VectorNormalize(edge02);

        float dot = x87f_store_f32(x87f_add(
            x87f_add(x87f_mul(x87f_load_f32(edge01[0]), x87f_load_f32(edge02[0])),
                     x87f_mul(x87f_load_f32(edge01[1]),
                              x87f_load_f32(edge02[1]))),
            x87f_mul(x87f_load_f32(edge01[2]), x87f_load_f32(edge02[2]))));
        /* inverseDot = 1 / (1 - dot*dot) */
        float inverseDot = x87f_store_f32(x87f_div(
            x87f_load_f32(1.0f),
            x87f_sub(x87f_load_f32(1.0f),
                     x87f_mul(x87f_load_f32(dot), x87f_load_f32(dot)))));

        for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
            edge01[axis] = x87f_store_f32(
                x87f_mul(x87f_load_f32(edge01[axis]), x87f_load_f32(inverseDot)));
            edge02[axis] = x87f_store_f32(
                x87f_mul(x87f_load_f32(edge02[axis]), x87f_load_f32(inverseDot)));
        }

        /* Each svec/tvec component is TWO stores: the numerator
         * edgeA[k] + (-dot)*edgeB[k] is rounded to a float slot (fchs;fmul;
         * fld;faddp;fstp), then that reloaded float is multiplied by 1/length
         * (fld1;fdiv length; fmulp; fstp) — the 1/length is formed inline per
         * component. */
        for (int32_t k = 0; k < CM_PATCH_VECTOR_AXIS_COUNT; ++k) {
            float numerator = x87f_store_f32(
                x87f_add(x87f_mul(x87f_neg(x87f_load_f32(dot)),
                                  x87f_load_f32(edge02[k])),
                         x87f_load_f32(edge01[k])));
            triangle->svec[k] = x87f_store_f32(x87f_mul(
                x87f_load_f32(numerator),
                x87f_div(x87f_load_f32(1.0f), x87f_load_f32(edge01Length))));
        }
        triangle->svec[3] = x87f_store_f32(x87f_add(
            x87f_add(x87f_mul(x87f_load_f32(triangle->svec[0]),
                              x87f_load_f32(vertices[first][0])),
                     x87f_mul(x87f_load_f32(triangle->svec[1]),
                              x87f_load_f32(vertices[first][1]))),
            x87f_mul(x87f_load_f32(triangle->svec[2]),
                     x87f_load_f32(vertices[first][2]))));

        for (int32_t k = 0; k < CM_PATCH_VECTOR_AXIS_COUNT; ++k) {
            float numerator = x87f_store_f32(
                x87f_add(x87f_mul(x87f_neg(x87f_load_f32(dot)),
                                  x87f_load_f32(edge01[k])),
                         x87f_load_f32(edge02[k])));
            triangle->tvec[k] = x87f_store_f32(x87f_mul(
                x87f_load_f32(numerator),
                x87f_div(x87f_load_f32(1.0f), x87f_load_f32(edge02Length))));
        }
        triangle->tvec[3] = x87f_store_f32(x87f_add(
            x87f_add(x87f_mul(x87f_load_f32(triangle->tvec[0]),
                              x87f_load_f32(vertices[first][0])),
                     x87f_mul(x87f_load_f32(triangle->tvec[1]),
                              x87f_load_f32(vertices[first][1]))),
            x87f_mul(x87f_load_f32(triangle->tvec[2]),
                     x87f_load_f32(vertices[first][2]))));
#else
        vec3_t edge01 = {
            vertices[second][0] - vertices[first][0],
            vertices[second][1] - vertices[first][1],
            vertices[second][2] - vertices[first][2],
        };
        float edge01Length = (float)VectorNormalize(edge01);
        vec3_t edge02 = {
            vertices[third][0] - vertices[first][0],
            vertices[third][1] - vertices[first][1],
            vertices[third][2] - vertices[first][2],
        };
        float edge02Length = (float)VectorNormalize(edge02);

        float dot = (edge01[0] * edge02[0]) + (edge01[1] * edge02[1]) +
                    (edge01[2] * edge02[2]);
        float inverseDot = 1.0f / (1.0f - (dot * dot));

        for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
            edge01[axis] *= inverseDot;
            edge02[axis] *= inverseDot;
        }

        /* The original rounds the (edgeA - dot*edgeB) numerator to a float
         * before multiplying by 1/length (it spills it to the svec/tvec slot
         * and reloads). A natural `(edgeA - dot*edgeB) * (1/length)` keeps the
         * numerator in an 80-bit register on both x87 and -mfpmath=387 at -O0
         * and does NOT match the stock binary, so use explicit float temps. */
        for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
            float numerator = edge01[axis] - (dot * edge02[axis]);
            triangle->svec[axis] = numerator * (1.0f / edge01Length);
        }
        triangle->svec[3] =
            (triangle->svec[0] * vertices[first][0]) +
            (triangle->svec[1] * vertices[first][1]) +
            (triangle->svec[2] * vertices[first][2]);

        for (int32_t axis = 0; axis < CM_PATCH_VECTOR_AXIS_COUNT; ++axis) {
            float numerator = edge02[axis] - (dot * edge01[axis]);
            triangle->tvec[axis] = numerator * (1.0f / edge02Length);
        }
        triangle->tvec[3] =
            (triangle->tvec[0] * vertices[first][0]) +
            (triangle->tvec[1] * vertices[first][1]) +
            (triangle->tvec[2] * vertices[first][2]);
#endif

        int16_t remap = cm_modelVertexRemap[first];
        triangle->vertices[0] =
            remap < 0 ? NULL : &collisionVertices[remap];
        remap = cm_modelVertexRemap[second];
        triangle->vertices[1] =
            remap < 0 ? NULL : &collisionVertices[remap];
        remap = cm_modelVertexRemap[third];
        triangle->vertices[2] =
            remap < 0 ? NULL : &collisionVertices[remap];

        int32_t edgeIndex =
            FindEdge((int32_t)edgeCount, cm_modelEdges, second, third);
        triangle->oppositeEdges[0] =
            edgeIndex < 0
                ? NULL
                : &collisionEdges[edgeIndex];
        edgeIndex = FindEdge((int32_t)edgeCount, cm_modelEdges, third, first);
        triangle->oppositeEdges[1] =
            edgeIndex < 0
                ? NULL
                : &collisionEdges[edgeIndex];
        edgeIndex = FindEdge((int32_t)edgeCount, cm_modelEdges, first, second);
        triangle->oppositeEdges[2] =
            edgeIndex < 0
                ? NULL
                : &collisionEdges[edgeIndex];
    }

    soup->negativeZOnly =
        hasNegativeZPlane != qfalse && hasPositiveZPlane == qfalse ? qtrue
                                                                   : qfalse;

    ClearBounds(bounds[0], bounds[1]);
    for (uint32_t vertexIndex = 0; vertexIndex < vertexCount; ++vertexIndex) {
        AddPointToBounds(vertices[vertexIndex], bounds[0], bounds[1]);

        int16_t remap = cm_modelVertexRemap[vertexIndex];
        if (remap >= 0) {
            collisionVertices[remap].checkcount = 0;
            collisionVertices[remap].position[0] = vertices[vertexIndex][0];
            collisionVertices[remap].position[1] = vertices[vertexIndex][1];
            collisionVertices[remap].position[2] = vertices[vertexIndex][2];
        }
    }

    for (int32_t edgeIndex = 0; edgeIndex < (int32_t)edgeCount; ++edgeIndex) {
        coduo_cm_collision_edge_t *edge = &collisionEdges[edgeIndex];
        int16_t first = cm_modelEdges[edgeIndex].first;
        int16_t second = cm_modelEdges[edgeIndex].second;

        edge->checkcount = 0;
        edge->origin[0] = vertices[first][0];
        edge->origin[1] = vertices[first][1];
        edge->origin[2] = vertices[first][2];
        edge->unitDirection[0] = vertices[second][0] - vertices[first][0];
        edge->unitDirection[1] = vertices[second][1] - vertices[first][1];
        edge->unitDirection[2] = vertices[second][2] - vertices[first][2];
        edge->length = (float)VectorNormalize(edge->unitDirection);
        PerpendicularVector(edge->radialAxes[0], edge->unitDirection);
        CrossProduct(edge->unitDirection, edge->radialAxes[0],
                     edge->radialAxes[1]);
    }

    return (coduo_patch_collide_t *)soup;
}
