1#include <math.h>
2#include <stdint.h>
3
4double remquo(double x, double y, int* quo) {
5    union {
6        double f;
7        uint64_t i;
8    } ux = {x}, uy = {y};
9    int ex = ux.i >> 52 & 0x7ff;
10    int ey = uy.i >> 52 & 0x7ff;
11    int sx = ux.i >> 63;
12    int sy = uy.i >> 63;
13    uint32_t q;
14    uint64_t i;
15    uint64_t uxi = ux.i;
16
17    *quo = 0;
18    if (uy.i << 1 == 0 || isnan(y) || ex == 0x7ff)
19        return (x * y) / (x * y);
20    if (ux.i << 1 == 0)
21        return x;
22
23    /* normalize x and y */
24    if (!ex) {
25        for (i = uxi << 12; i >> 63 == 0; ex--, i <<= 1)
26            ;
27        uxi <<= -ex + 1;
28    } else {
29        uxi &= -1ULL >> 12;
30        uxi |= 1ULL << 52;
31    }
32    if (!ey) {
33        for (i = uy.i << 12; i >> 63 == 0; ey--, i <<= 1)
34            ;
35        uy.i <<= -ey + 1;
36    } else {
37        uy.i &= -1ULL >> 12;
38        uy.i |= 1ULL << 52;
39    }
40
41    q = 0;
42    if (ex < ey) {
43        if (ex + 1 == ey)
44            goto end;
45        return x;
46    }
47
48    /* x mod y */
49    for (; ex > ey; ex--) {
50        i = uxi - uy.i;
51        if (i >> 63 == 0) {
52            uxi = i;
53            q++;
54        }
55        uxi <<= 1;
56        q <<= 1;
57    }
58    i = uxi - uy.i;
59    if (i >> 63 == 0) {
60        uxi = i;
61        q++;
62    }
63    if (uxi == 0)
64        ex = -60;
65    else
66        for (; uxi >> 52 == 0; uxi <<= 1, ex--)
67            ;
68end:
69    /* scale result and decide between |x| and |x|-|y| */
70    if (ex > 0) {
71        uxi -= 1ULL << 52;
72        uxi |= (uint64_t)ex << 52;
73    } else {
74        uxi >>= -ex + 1;
75    }
76    ux.i = uxi;
77    x = ux.f;
78    if (sy)
79        y = -y;
80    if (ex == ey || (ex + 1 == ey && (2 * x > y || (2 * x == y && q % 2)))) {
81        x -= y;
82        q++;
83    }
84    q &= 0x7fffffff;
85    *quo = sx ^ sy ? -(int)q : (int)q;
86    return sx ? -x : x;
87}
88