1#include "libm.h"
2
3#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
4long double fmodl(long double x, long double y) {
5    return fmod(x, y);
6}
7#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
8long double fmodl(long double x, long double y) {
9    union ldshape ux = {x}, uy = {y};
10    int ex = ux.i.se & 0x7fff;
11    int ey = uy.i.se & 0x7fff;
12    int sx = ux.i.se & 0x8000;
13
14    if (y == 0 || isnan(y) || ex == 0x7fff)
15        return (x * y) / (x * y);
16    ux.i.se = ex;
17    uy.i.se = ey;
18    if (ux.f <= uy.f) {
19        if (ux.f == uy.f)
20            return 0 * x;
21        return x;
22    }
23
24    /* normalize x and y */
25    if (!ex) {
26        ux.f *= 0x1p120f;
27        ex = ux.i.se - 120;
28    }
29    if (!ey) {
30        uy.f *= 0x1p120f;
31        ey = uy.i.se - 120;
32    }
33
34/* x mod y */
35#if LDBL_MANT_DIG == 64
36    uint64_t i, zx, my;
37    zx = ux.i.m;
38    my = uy.i.m;
39    for (; ex > ey; ex--) {
40        i = zx - my;
41        if (zx >= my) {
42            if (i == 0)
43                return 0 * x;
44            zx = 2 * i;
45        } else if (2 * zx < zx) {
46            zx = 2 * zx - my;
47        } else {
48            zx = 2 * zx;
49        }
50    }
51    i = zx - my;
52    if (zx >= my) {
53        if (i == 0)
54            return 0 * x;
55        zx = i;
56    }
57    for (; zx >> 63 == 0; zx *= 2, ex--)
58        ;
59    ux.i.m = zx;
60#elif LDBL_MANT_DIG == 113
61    uint64_t hi, lo, xhi, xlo, yhi, ylo;
62    xhi = (ux.i2.hi & -1ULL >> 16) | 1ULL << 48;
63    yhi = (uy.i2.hi & -1ULL >> 16) | 1ULL << 48;
64    xlo = ux.i2.lo;
65    ylo = uy.i2.lo;
66    for (; ex > ey; ex--) {
67        hi = xhi - yhi;
68        lo = xlo - ylo;
69        if (xlo < ylo)
70            hi -= 1;
71        if (hi >> 63 == 0) {
72            if ((hi | lo) == 0)
73                return 0 * x;
74            xhi = 2 * hi + (lo >> 63);
75            xlo = 2 * lo;
76        } else {
77            xhi = 2 * xhi + (xlo >> 63);
78            xlo = 2 * xlo;
79        }
80    }
81    hi = xhi - yhi;
82    lo = xlo - ylo;
83    if (xlo < ylo)
84        hi -= 1;
85    if (hi >> 63 == 0) {
86        if ((hi | lo) == 0)
87            return 0 * x;
88        xhi = hi;
89        xlo = lo;
90    }
91    for (; xhi >> 48 == 0; xhi = 2 * xhi + (xlo >> 63), xlo = 2 * xlo, ex--)
92        ;
93    ux.i2.hi = xhi;
94    ux.i2.lo = xlo;
95#endif
96
97    /* scale result */
98    if (ex <= 0) {
99        ux.i.se = (ex + 120) | sx;
100        ux.f *= 0x1p-120f;
101    } else
102        ux.i.se = ex | sx;
103    return ux.f;
104}
105#endif
106