#include <math.h>

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

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

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

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

    if (startDistanceMinusRadius <= 0.0f) {
        traceWork->trace.fraction = 0.0f;
        traceWork->trace.startsolid = 1;
        VectorNormalize2(delta, traceWork->trace.normal);
        traceWork->trace.contents = CM_TempBoxModelContents();

        delta[0] = end[0] - sphereOrigin[0];
        delta[1] = end[1] - sphereOrigin[1];
        delta[2] = end[2] - sphereOrigin[2];
        /* delta.delta 3-mul/2-add dot kept 80-bit, full-width compare vs
         * radiusSquared (inline, not stored) -> shim. */
#if EMULATE_X87
        if (x87f_le(
                x87f_add(x87f_add(
                    x87f_mul(x87f_load_f32(delta[0]), x87f_load_f32(delta[0])),
                    x87f_mul(x87f_load_f32(delta[1]), x87f_load_f32(delta[1]))),
                    x87f_mul(x87f_load_f32(delta[2]), x87f_load_f32(delta[2]))),
                x87f_load_f32(radiusSquared))) {
            traceWork->trace.allsolid = 1;
        }
#else
        if ((((delta[0] * delta[0]) + (delta[1] * delta[1])) +
             (delta[2] * delta[2])) <= radiusSquared) {
            traceWork->trace.allsolid = 1;
        }
#endif

        return 0;
    }

    /* traceWork->delta . delta dot kept 80-bit, stored to float -> shim. */
#if EMULATE_X87
    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]))));
#else
    const float deltaDot =
        ((traceWork->delta[0] * delta[0]) +
         (traceWork->delta[1] * delta[1])) +
        (traceWork->delta[2] * delta[2]);
#endif

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

    const float deltaLengthSquared = traceWork->deltaLengthSquared;
    /* deltaDot^2 - lenSq*startDist kept 80-bit, stored to float (0x805631d). */
#if EMULATE_X87
    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))));
#else
    const float discriminant =
        (deltaDot * deltaDot) -
        (deltaLengthSquared * startDistanceMinusRadius);
#endif

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

    const float startDistance = VectorNormalize2(delta, normal);
    /* (-deltaDot - sqrt((double)disc)) / lenSq kept 80-bit (0x8056349: fchs;
     * fsub deltaDot; fdiv lenSq); enterFraction = base + (startDist*0.125)/
     * deltaDot, both one store -> shim. */
#if EMULATE_X87
    const float baseFraction = x87f_store_f32(x87f_div(
        x87f_sub(x87f_neg(x87f_load_f32(deltaDot)),
                 x87f_load_f64(sqrt((double)discriminant))),
        x87f_load_f32(deltaLengthSquared)));
    const float enterFraction = x87f_store_f32(x87f_add(
        x87f_load_f32(baseFraction),
        x87f_div(x87f_mul(x87f_load_f32(startDistance), x87f_load_f32(0.125f)),
                 x87f_load_f32(deltaDot))));
#else
    const float baseFraction =
        (-deltaDot - sqrt((double)discriminant)) / deltaLengthSquared;
    const float enterFraction =
        baseFraction + ((startDistance * 0.125f) / deltaDot);
#endif

    if (enterFraction < traceWork->trace.fraction) {
        if (enterFraction > 0.0f) {
            traceWork->trace.fraction = enterFraction;
        } else {
            traceWork->trace.fraction = 0.0f;
        }

        traceWork->trace.normal[0] = normal[0];
        traceWork->trace.normal[1] = normal[1];
        traceWork->trace.normal[2] = normal[2];
        traceWork->trace.contents = CM_TempBoxModelContents();

        return 0;
    }

    return 1;
}
