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

#include <math.h>

qboolean CM_SightTraceSphereThroughSphere(
    const traceWork_t *traceWork,
    const vec3_t start,
    const vec3_t end,
    const vec3_t sphereOrigin,
    float sphereRadius)
{
    vec3_t delta;
    vec3_t normal;

    (void)end;

    delta[0] = start[0] - sphereOrigin[0];
    delta[1] = start[1] - sphereOrigin[1];
    delta[2] = start[2] - sphereOrigin[2];

    const float expandedRadiusSquared =
        (sphereRadius + traceWork->sphere.radius) *
        (sphereRadius + traceWork->sphere.radius);
    const float startDistanceMinusRadius =
        ((delta[0] * delta[0] + delta[1] * delta[1]) +
         delta[2] * delta[2]) -
        expandedRadiusSquared;

    if (startDistanceMinusRadius <= 0.0f) {
        return 0;
    }

#if EMULATE_X87
    /* 3-comp dot, discriminant, rootFraction, enterFraction each 80-bit one
     * store -> shim. */
    const float deltaDot = x87f_store_f32(x87f_add(
        x87f_add(x87f_mul(x87f_load_f32(traceWork->delta[0]), x87f_load_f32(delta[0])),
                 x87f_mul(x87f_load_f32(traceWork->delta[1]), x87f_load_f32(delta[1]))),
        x87f_mul(x87f_load_f32(traceWork->delta[2]), x87f_load_f32(delta[2]))));

    if (0.0f <= deltaDot) {
        return 1;
    }

    const float deltaLengthSquared = traceWork->deltaLengthSquared;
    const float discriminant = x87f_store_f32(x87f_sub(
        x87f_mul(x87f_load_f32(deltaDot), x87f_load_f32(deltaDot)),
        x87f_mul(x87f_load_f32(deltaLengthSquared),
                 x87f_load_f32(startDistanceMinusRadius))));

    if (discriminant < 0.0f) {
        return 1;
    }

    const float startDistance = VectorNormalize2(delta, normal);
    const float rootFraction = x87f_store_f32(x87f_div(
        x87f_sub(x87f_neg(x87f_load_f32(deltaDot)),
                 x87f_load_f64(sqrt((double)discriminant))),
        x87f_load_f32(deltaLengthSquared)));
    /*
     * ORIGINAL_BINARY_BUG: stock multiplies deltaDot by the clip epsilon and
     * divides by startDistance, the reciprocal of the dimensionless
     * sphere-trace adjustment.  See
     * original-bugs/coduo-sight-capsule-reciprocal-epsilon.md.
     */
    const float enterFraction = x87f_store_f32(x87f_add(
        x87f_load_f32(rootFraction),
        x87f_div(x87f_mul(x87f_load_f32(deltaDot), x87f_load_f32(0.125f)),
                 x87f_load_f32(startDistance))));
#else
    const float deltaDot =
        (traceWork->delta[0] * delta[0] +
         traceWork->delta[1] * delta[1]) +
        traceWork->delta[2] * delta[2];

    if (0.0f <= deltaDot) {
        return 1;
    }

    const float deltaLengthSquared = traceWork->deltaLengthSquared;
    const float discriminant =
        deltaDot * deltaDot -
        deltaLengthSquared * startDistanceMinusRadius;

    if (discriminant < 0.0f) {
        return 1;
    }

    const float startDistance = VectorNormalize2(delta, normal);
    const float rootFraction =
        (-deltaDot - sqrt((double)discriminant)) / deltaLengthSquared;
    /* ORIGINAL_BINARY_BUG: preserved reciprocal epsilon adjustment above. */
    const float enterFraction =
        rootFraction + (deltaDot * 0.125f) / startDistance;
#endif

    if (enterFraction < traceWork->trace.fraction) {
        return 0;
    }

    return 1;
}
