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