#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <math.h>

#ifdef HAVE_CONFIG_H
#include <errno.h>
#endif

/*
 * The floating-point remainder x - n*y with n = x/y truncated toward zero
 * (C23 7.12.10.1), all three precisions.  The remainder is always exact
 * and carries x's sign (a zero result is a signed zero with x's sign), so
 * every implementation must agree bit for bit; the result is computed with
 * a restoring long division on the significands, in the style of the musl
 * fmod family but derived here from first principles.
 *
 * Each operand is normalized into
 *
 *     value = m * 2^(ex - W)
 *
 * where m is the significand with its msb pinned to bit P (P = 52/23/63
 * and W = 1075/150/16446 for double/float/80-bit) and ex is the exponent
 * in the same units as the biased field (so ex is the stored exponent for
 * normals and goes negative for subnormals).  Aligned this way, one binary
 * long-division step per exponent difference decides whether a multiple of
 * the divisor fits: the current remainder significand is compared against
 * the divisor significand, the divisor is subtracted once when it fits,
 * and the remainder is doubled for the next, half-weight, digit.  The
 * comparison invariant keeps the remainder below twice the divisor at
 * every step, so one subtraction per digit is always enough and every
 * subtraction (and the exact zero test) is a plain integer operation.
 *
 * Domain errors mirror the host glibc: fmod(+-0, +-0), fmod(+-Inf, y) and
 * any NaN argument return a NaN (the classic indefinite pattern, which is
 * what glibc's x87 path produces) and, in the library build only, set
 * errno to EDOM.  fmod(x, +-Inf) is x for a finite x, and |x| <= |y| hands
 * x back unchanged, which keeps the signed-zero and exact cases exact.
 */

/* The double format: bit 63 the sign, bits 62..52 the exponent biased by
 * 1023, bits 51..0 the fraction.  Value = m * 2^(ex - 1075) with the msb
 * of m at bit 52. */
static double
fmod_d(double x, double y) // NOLINT(bugprone-easily-swappable-parameters)
{
    unsigned long long xb;
    unsigned long long yb;
    unsigned long long sx;
    unsigned long long mx;
    unsigned long long my;
    unsigned long long s;
    unsigned long long i;
    int ex;
    int ey;

    __builtin_memcpy(&xb, &x, sizeof xb);
    __builtin_memcpy(&yb, &y, sizeof yb);
    sx = xb & (1ULL << 63);
    ex = (int)((xb >> 52) & 0x7ff);
    ey = (int)((yb >> 52) & 0x7ff);

    /* x not finite, y zero, or y a NaN: domain error. */
    if (ex == 0x7ff || (yb << 1) == 0 || (ey == 0x7ff && (yb & 0xFFFFFFFFFFFFFULL) != 0))
    {
#ifdef HAVE_CONFIG_H
        errno = EDOM;
#endif
        return __builtin_nan(""); /* indefinite NaN, sign negative */
    }
    /* x zero or a finite x below an infinite y: x unchanged.  (x a NaN or
     * +-Inf fell into the domain branch above, which returns a NaN for
     * both.) */
    if ((xb << 1) == 0 || ey == 0x7ff)
    {
        return x;
    }

    /* Compare magnitudes directly on the (biased exponent, fraction)
     * words: finite IEEE magnitudes order lexicographically, so a plain
     * integer compare of the sign-stripped words decides |x| vs |y|. */
    if ((xb & 0x7FFFFFFFFFFFFFFFULL) <= (yb & 0x7FFFFFFFFFFFFFFFULL))
    {
        if ((xb & 0x7FFFFFFFFFFFFFFFULL) == (yb & 0x7FFFFFFFFFFFFFFFULL))
        {
            return x * 0.0; /* exact: result is +-0 with x's sign */
        }
        return x;
    }

    /* Normalize x: subnormals (ex == 0) get their msb shifted up to bit
     * 52 and ex counts the shift below the smallest normal. */
    mx = xb & 0xFFFFFFFFFFFFFULL;
    if (ex == 0)
    {
        s = 63 - (unsigned long long)__builtin_clzll(mx);
        mx <<= (52 - s);
        ex = (int)s - 51;
    }
    else
    {
        mx |= 1ULL << 52;
    }
    my = yb & 0xFFFFFFFFFFFFFULL;
    if (ey == 0)
    {
        s = 63 - (unsigned long long)__builtin_clzll(my);
        my <<= (52 - s);
        ey = (int)s - 51;
    }
    else
    {
        my |= 1ULL << 52;
    }

    /* Long division: subtract the aligned divisor significand once per
     * bit of quotient, doubling the remainder between bits. */
    for (; ex > ey; ex--)
    {
        if (mx >= my)
        {
            i = mx - my;
            if (i == 0)
            {
                return x * 0.0;
            }
            mx = i;
        }
        mx <<= 1;
    }
    i = mx - my;
    if (mx >= my)
    {
        if (i == 0)
        {
            return x * 0.0;
        }
        mx = i;
    }

    /* Bring the msb of the remainder back to bit 52. */
    for (; (mx >> 52) == 0; mx <<= 1, ex--)
    {
    }

    if (ex > 0)
    {
        xb = sx | ((unsigned long long)ex << 52) | (mx - (1ULL << 52));
    }
    else
    {
        xb = sx | (mx >> (1 - ex));
    }
    __builtin_memcpy(&x, &xb, sizeof xb);
    return x;
}

/* The float format: bit 31 the sign, bits 30..23 the exponent biased by
 * 127, bits 22..0 the fraction.  Value = m * 2^(ex - 150) with the msb of
 * m at bit 23. */
static float
fmod_f(float x, float y) // NOLINT(bugprone-easily-swappable-parameters)
{
    unsigned int xb;
    unsigned int yb;
    unsigned int sx;
    unsigned int mx;
    unsigned int my;
    unsigned int s;
    unsigned int i;
    int ex;
    int ey;

    __builtin_memcpy(&xb, &x, sizeof xb);
    __builtin_memcpy(&yb, &y, sizeof yb);
    sx = xb & (1U << 31);
    ex = (int)((xb >> 23) & 0xff);
    ey = (int)((yb >> 23) & 0xff);

    if (ex == 0xff || (yb << 1) == 0 || (ey == 0xff && (yb & 0x7FFFFFU) != 0))
    {
#ifdef HAVE_CONFIG_H
        errno = EDOM;
#endif
        return __builtin_nanf("");
    }
    if ((xb << 1) == 0 || ey == 0xff)
    {
        return x;
    }

    if ((xb & 0x7FFFFFFFU) <= (yb & 0x7FFFFFFFU))
    {
        if ((xb & 0x7FFFFFFFU) == (yb & 0x7FFFFFFFU))
        {
            return x * 0.0f;
        }
        return x;
    }

    mx = xb & 0x7FFFFFU;
    if (ex == 0)
    {
        s = 31 - (unsigned int)__builtin_clz(mx);
        mx <<= (23 - s);
        ex = (int)s - 22;
    }
    else
    {
        mx |= 1U << 23;
    }
    my = yb & 0x7FFFFFU;
    if (ey == 0)
    {
        s = 31 - (unsigned int)__builtin_clz(my);
        my <<= (23 - s);
        ey = (int)s - 22;
    }
    else
    {
        my |= 1U << 23;
    }

    for (; ex > ey; ex--)
    {
        if (mx >= my)
        {
            i = mx - my;
            if (i == 0)
            {
                return x * 0.0f;
            }
            mx = i;
        }
        mx <<= 1;
    }
    i = mx - my;
    if (mx >= my)
    {
        if (i == 0)
        {
            return x * 0.0f;
        }
        mx = i;
    }

    for (; (mx >> 23) == 0; mx <<= 1, ex--)
    {
    }

    if (ex > 0)
    {
        xb = sx | ((unsigned int)ex << 23) | (mx - (1U << 23));
    }
    else
    {
        xb = sx | (mx >> (1 - ex));
    }
    __builtin_memcpy(&x, &xb, sizeof xb);
    return x;
}

/*
 * The x86 80-bit extended format: 64-bit significand m with an explicit
 * integer bit in bytes 0..7, sign/exponent word se in bytes 8..9.  Value =
 * m * 2^(ex - 16446) with the msb of m at bit 63; a canonical normal has
 * ex = se & 0x7fff, and a subnormal (se field 0) or a defensive unnormal
 * (nonzero field with m < 2^63) is normalized by shifting m's msb up to
 * bit 63.  The division reuses the same digit loop as the narrower
 * formats; only the doubling step differs, because a full-width
 * significand has no free bit above its msb (see the three-way branch in
 * the loop below).
 */
struct fmod_ld_word
{
    unsigned long long m;
    unsigned short se;
};

static long double
fmod_l(long double x, long double y) // NOLINT(bugprone-easily-swappable-parameters)
{
    struct fmod_ld_word p;
    struct fmod_ld_word q;
    unsigned long long s;
    unsigned long long mx;
    unsigned long long my;
    unsigned long long i;
    int ex;
    int ey;

    __builtin_memcpy(&p, &x, sizeof p);
    __builtin_memcpy(&q, &y, sizeof q);
    ex = p.se & 0x7fff;
    ey = q.se & 0x7fff;

    /* x not finite, y zero, or y a NaN: domain error. */
    if (ex == 0x7fff || (q.m == 0 && ey == 0) ||
        (ey == 0x7fff && (q.m & 0x7FFFFFFFFFFFFFFFULL) != 0))
    {
#ifdef HAVE_CONFIG_H
        errno = EDOM;
#endif
        p.m = 0xC000000000000000ULL;
        p.se = 0xFFFF;
        __builtin_memcpy(&x, &p, sizeof p);
        return x;
    }
    /* x zero or a finite x below an infinite y: x unchanged.  (x a NaN or
     * +-Inf fell into the domain branch above.) */
    if ((p.m == 0 && ex == 0) || ey == 0x7fff)
    {
        return x;
    }

    /* Normalize both significands to the [2^63, 2^64) frame; a subnormal
     * (field 0) or defensive unnormal (field nonzero, msb below 63) has
     * its msb shifted up to bit 63 with the exponent adjusted, so every
     * value below obeys value = m * 2^(ex - 16446). */
    if (p.m < 0x8000000000000000ULL)
    {
        s = 63 - (unsigned long long)__builtin_clzll(p.m);
        p.m <<= (63 - s);
        ex = (ex == 0 ? 1 : ex) + (int)s - 63;
    }
    if (q.m < 0x8000000000000000ULL)
    {
        s = 63 - (unsigned long long)__builtin_clzll(q.m);
        q.m <<= (63 - s);
        ey = (ey == 0 ? 1 : ey) + (int)s - 63;
    }
    mx = p.m;
    my = q.m;

    if (ex < ey || (ex == ey && mx <= my))
    {
        if (ex == ey && mx == my)
        {
            return x * 0.0L;
        }
        return x;
    }

    /* Long division with a full-width significand: a 64-bit divisor
     * leaves no headroom above its own msb for the per-bit doubling, so
     * the doubled remainder that overflows is exactly one divisor at the
     * next, half, weight and is absorbed by a subtraction there. */
    for (; ex > ey; ex--)
    {
        i = mx - my;
        if (mx >= my)
        {
            if (i == 0)
            {
                return x * 0.0L;
            }
            mx = 2 * i;
        }
        else if (2 * mx < mx)
        {
            mx = 2 * mx - my;
        }
        else
        {
            mx = 2 * mx;
        }
    }
    i = mx - my;
    if (mx >= my)
    {
        if (i == 0)
        {
            return x * 0.0L;
        }
        mx = i;
    }

    /* Bring the msb of the remainder back to bit 63. */
    for (; mx < 0x8000000000000000ULL; mx <<= 1, ex--)
    {
    }

    if (ex > 0)
    {
        p.se = (unsigned short)((p.se & 0x8000) | (unsigned short)ex);
        p.m = mx;
    }
    else
    {
        p.se = (unsigned short)(p.se & 0x8000);
        p.m = mx >> (1 - ex);
    }
    __builtin_memcpy(&x, &p, sizeof p);
    return x;
}

/*
 * The float remainder, as fmod_d.
 */
float
fmodf(float x, float y) // NOLINT(bugprone-easily-swappable-parameters)
{
    return fmod_f(x, y);
}

/*
 * The double remainder, as fmod_d.
 */
double
fmod(double x, double y) // NOLINT(bugprone-easily-swappable-parameters)
{
    return fmod_d(x, y);
}

/*
 * The long double remainder, as fmod_d.
 */
long double
fmodl(long double x, long double y) // NOLINT(bugprone-easily-swappable-parameters)
{
    return fmod_l(x, y);
}