1#include <math.h>
2#include <stdint.h>
3
4float remquof(float x, float y, int* quo) {
5    union {
6        float f;
7        uint32_t i;
8    } ux = {x}, uy = {y};
9    int ex = ux.i >> 23 & 0xff;
10    int ey = uy.i >> 23 & 0xff;
11    int sx = ux.i >> 31;
12    int sy = uy.i >> 31;
13    uint32_t q;
14    uint32_t i;
15    uint32_t uxi = ux.i;
16
17    *quo = 0;
18    if (uy.i << 1 == 0 || isnan(y) || ex == 0xff)
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 << 9; i >> 31 == 0; ex--, i <<= 1)
26            ;
27        uxi <<= -ex + 1;
28    } else {
29        uxi &= -1U >> 9;
30        uxi |= 1U << 23;
31    }
32    if (!ey) {
33        for (i = uy.i << 9; i >> 31 == 0; ey--, i <<= 1)
34            ;
35        uy.i <<= -ey + 1;
36    } else {
37        uy.i &= -1U >> 9;
38        uy.i |= 1U << 23;
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 >> 31 == 0) {
52            uxi = i;
53            q++;
54        }
55        uxi <<= 1;
56        q <<= 1;
57    }
58    i = uxi - uy.i;
59    if (i >> 31 == 0) {
60        uxi = i;
61        q++;
62    }
63    if (uxi == 0)
64        ex = -30;
65    else
66        for (; uxi >> 23 == 0; uxi <<= 1, ex--)
67            ;
68end:
69    /* scale result and decide between |x| and |x|-|y| */
70    if (ex > 0) {
71        uxi -= 1U << 23;
72        uxi |= (uint32_t)ex << 23;
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