#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double fmodl(long double x, long double y)
{
    return fmod(x, y);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
long double fmodl(long double x, long double y)
{
    union ldshape ux = {x}, uy = {y};
    int ex = ux.i.se & 0x7fff;
    int ey = uy.i.se & 0x7fff;
    int sx = ux.i.se & 0x8000;

    if (y == 0 || isnan(y) || ex == 0x7fff)
        return (x * y) / (x * y);
    ux.i.se = ex;
    uy.i.se = ey;
    if (ux.f <= uy.f)
    {
        if (ux.f == uy.f)
            return 0 * x;
        return x;
    }

    /* normalize x and y */
    if (!ex)
    {
        ux.f *= 0x1p120f;
        ex = ux.i.se - 120;
    }
    if (!ey)
    {
        uy.f *= 0x1p120f;
        ey = uy.i.se - 120;
    }

    /* x mod y */
#    if LDBL_MANT_DIG == 64
    uint64_t i, mx, my;
    mx = ux.i.m;
    my = uy.i.m;
    for (; ex > ey; ex--)
    {
        i = mx - my;
        if (mx >= my)
        {
            if (i == 0)
                return 0 * x;
            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 0 * x;
        mx = i;
    }
    for (; mx >> 63 == 0; mx *= 2, ex--)
        ;
    ux.i.m = mx;
#    elif LDBL_MANT_DIG == 113
    uint64_t hi, lo, xhi, xlo, yhi, ylo;
    xhi = (ux.i2.hi & -1ULL >> 16) | 1ULL << 48;
    yhi = (uy.i2.hi & -1ULL >> 16) | 1ULL << 48;
    xlo = ux.i2.lo;
    ylo = uy.i2.lo;
    for (; ex > ey; ex--)
    {
        hi = xhi - yhi;
        lo = xlo - ylo;
        if (xlo < ylo)
            hi -= 1;
        if (hi >> 63 == 0)
        {
            if ((hi | lo) == 0)
                return 0 * x;
            xhi = 2 * hi + (lo >> 63);
            xlo = 2 * lo;
        }
        else
        {
            xhi = 2 * xhi + (xlo >> 63);
            xlo = 2 * xlo;
        }
    }
    hi = xhi - yhi;
    lo = xlo - ylo;
    if (xlo < ylo)
        hi -= 1;
    if (hi >> 63 == 0)
    {
        if ((hi | lo) == 0)
            return 0 * x;
        xhi = hi;
        xlo = lo;
    }
    for (; xhi >> 48 == 0; xhi = 2 * xhi + (xlo >> 63), xlo = 2 * xlo, ex--)
        ;
    ux.i2.hi = xhi;
    ux.i2.lo = xlo;
#    endif

    /* scale result */
    if (ex <= 0)
    {
        ux.i.se = (ex + 120) | sx;
        ux.f *= 0x1p-120f;
    }
    else
        ux.i.se = ex | sx;
    return ux.f;
}
#endif
