1#include <math.h>
2#include <stdint.h>
3
4float fmodf(float x, float y) {
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    uint32_t sx = ux.i & 0x80000000;
12    uint32_t i;
13    uint32_t uxi = ux.i;
14
15    if (uy.i << 1 == 0 || isnan(y) || ex == 0xff)
16        return (x * y) / (x * y);
17    if (uxi << 1 <= uy.i << 1) {
18        if (uxi << 1 == uy.i << 1)
19            return 0 * x;
20        return x;
21    }
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    /* x mod y */
42    for (; ex > ey; ex--) {
43        i = uxi - uy.i;
44        if (i >> 31 == 0) {
45            if (i == 0)
46                return 0 * x;
47            uxi = i;
48        }
49        uxi <<= 1;
50    }
51    i = uxi - uy.i;
52    if (i >> 31 == 0) {
53        if (i == 0)
54            return 0 * x;
55        uxi = i;
56    }
57    for (; uxi >> 23 == 0; uxi <<= 1, ex--)
58        ;
59
60    /* scale result up */
61    if (ex > 0) {
62        uxi -= 1U << 23;
63        uxi |= (uint32_t)ex << 23;
64    } else {
65        uxi >>= -ex + 1;
66    }
67    uxi |= sx;
68    ux.i = uxi;
69    return ux.f;
70}
71