#include <math.h>

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

qboolean CM_SightTraceCylinderThroughCylinder(
    const traceWork_t *traceWork,
    const vec3_t planeOrigin,
    float planeHalfHeight,
    float planeRadius)
{
    vec3_t startDelta;
    vec3_t hitNormal;

    startDelta[0] = traceWork->start[0] - planeOrigin[0];
    startDelta[1] = traceWork->start[1] - planeOrigin[1];
    startDelta[2] = traceWork->start[2] - planeOrigin[2];

    /* 0x805c0f8: the squared combined radius is rounded to its own float
     * slot before the radial dot subtracts it (as in the trace-side
     * cylinder probe). */
    const float combinedRadiusSquared =
        (planeRadius + traceWork->sphere.radius) *
        (planeRadius + traceWork->sphere.radius);
    const float radialStartDistance =
        (startDelta[0] * startDelta[0] + startDelta[1] * startDelta[1]) -
        combinedRadiusSquared;

    if (radialStartDistance <= 0.0f) {
        const float expandedHalfHeight =
            traceWork->sphere.halfheight - traceWork->sphere.radius +
            planeHalfHeight;

        if (expandedHalfHeight < startDelta[2] ||
            startDelta[2] < -expandedHalfHeight) {
            return 1;
        }

        return 0;
    }

#if EMULATE_X87
    /* 2-comp radial dot, discriminant, epsilon/base fractions each 80-bit one
     * store -> shim; hitFraction single add (native). */
    const float radialDeltaDot = x87f_store_f32(x87f_add(
        x87f_mul(x87f_load_f32(traceWork->delta[0]), x87f_load_f32(startDelta[0])),
        x87f_mul(x87f_load_f32(traceWork->delta[1]), x87f_load_f32(startDelta[1]))));
    if (0.0f <= radialDeltaDot) {
        return 1;
    }

    /*
     * ORIGINAL_BINARY_BUG: 0x805c1a8 loads the full XYZ delta length even
     * though this discriminant is the radial XY cylinder quadratic.  See
     * original-bugs/coduo-sight-cylinder-three-dimensional-quadratic.md.
     */
    const float deltaLengthSquared = traceWork->deltaLengthSquared;
    const float discriminant = x87f_store_f32(x87f_sub(
        x87f_mul(x87f_load_f32(radialDeltaDot), x87f_load_f32(radialDeltaDot)),
        x87f_mul(x87f_load_f32(deltaLengthSquared),
                 x87f_load_f32(radialStartDistance))));
    if (discriminant < 0.0f) {
        return 1;
    }

    startDelta[2] = 0.0f;
    const float radialStartLength = VectorNormalize2(startDelta, hitNormal);
    /*
     * ORIGINAL_BINARY_BUG: 0x805c1f8..0x805c206 computes
     * radialDeltaDot*0.125/radialStartLength, reciprocal to the dimensionless
     * trace-side adjustment.  See
     * original-bugs/coduo-sight-capsule-reciprocal-epsilon.md.
     */
    const float epsilonFraction = x87f_store_f32(x87f_div(
        x87f_mul(x87f_load_f32(radialDeltaDot), x87f_load_f32(0.125f)),
        x87f_load_f32(radialStartLength)));
    const float baseFraction = x87f_store_f32(x87f_div(
        x87f_sub(x87f_neg(x87f_load_f32(radialDeltaDot)),
                 x87f_load_f64(sqrt((double)discriminant))),
        x87f_load_f32(deltaLengthSquared)));
    const float hitFraction = baseFraction + epsilonFraction;
#else
    const float radialDeltaDot =
        traceWork->delta[0] * startDelta[0] +
        traceWork->delta[1] * startDelta[1];
    if (0.0f <= radialDeltaDot) {
        return 1;
    }

    const float deltaLengthSquared = traceWork->deltaLengthSquared;
    const float discriminant =
        radialDeltaDot * radialDeltaDot -
        deltaLengthSquared * radialStartDistance;
    if (discriminant < 0.0f) {
        return 1;
    }

    startDelta[2] = 0.0f;
    const float radialStartLength = VectorNormalize2(startDelta, hitNormal);
    /* ORIGINAL_BINARY_BUG: preserved reciprocal epsilon adjustment above. */
    const float epsilonFraction =
        (radialDeltaDot * 0.125f) / radialStartLength;
    const float baseFraction =
        (-radialDeltaDot - sqrt((double)discriminant)) / deltaLengthSquared;
    const float hitFraction = baseFraction + epsilonFraction;
#endif

    if (!(traceWork->trace.fraction > hitFraction)) {
        return 1;
    }

    const float expandedHalfHeight =
        traceWork->sphere.halfheight - traceWork->sphere.radius +
        planeHalfHeight;
    const float zAtHit =
        ((hitFraction - epsilonFraction) * traceWork->delta[2] +
         traceWork->start[2]) -
        planeOrigin[2];

    if (expandedHalfHeight < zAtHit || zAtHit < -expandedHalfHeight) {
        return 1;
    }

    return 0;
}
