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

#include "core_math_private.h"
#include "coduo_int32_bits.h"

long double Q_rsqrt(float value)
{
    uint32_t valueBits;
    uint32_t estimateBits;
    float estimate;

    memcpy(&valueBits, &value, sizeof(valueBits));
    estimateBits = UINT32_C(0x5f3759df) -
                   coduo_int32_sar_bits(valueBits, 1U);
    memcpy(&estimate, &estimateBits, sizeof(estimate));

#if EMULATE_X87
    /* x87-faithful transcription (stock coduo_lnxded 0x080660af): halfValue =
     * value*0.5f rounded to a float slot (fmulp; fstps). After the integer seed
     * (magic 0x5f3759df minus arithmetic-shifted bits), the single Newton step is
     * evaluated in 80-bit as y*(1.5f - ((halfValue*y)*y)) — the two y factors are
     * successive fmuls of the float slot, 1.5f applied via fsubrp, the final
     * product rounded to the float result (fstps). */
    const float halfValue =
        x87f_store_f32(x87f_mul(x87f_load_f32(value), x87f_load_f32(0.5f)));

    {
        const x87f y = x87f_load_f32(estimate);
        const x87f newton =
            x87f_sub(x87f_load_f32(1.5f),
                     x87f_mul(x87f_mul(x87f_load_f32(halfValue), y), y));
        estimate = x87f_store_f32(x87f_mul(y, newton));
    }
#else
    const float halfValue = value * 0.5f;

    estimate = estimate * (1.5f - halfValue * estimate * estimate);
#endif

    return estimate;
}
