1#include <math.h>
2#include <stdint.h>
3
4double fmod(double x, double y) {
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    uint64_t i;
13
14    /* in the followings uxi should be ux.i, but then gcc wrongly adds */
15    /* float load/store to inner loops ruining performance and code size */
16    uint64_t uxi = ux.i;
17
18    if (uy.i << 1 == 0 || isnan(y) || ex == 0x7ff)
19        return (x * y) / (x * y);
20    if (uxi << 1 <= uy.i << 1) {
21        if (uxi << 1 == uy.i << 1)
22            return 0 * x;
23        return x;
24    }
25
26    /* normalize x and y */
27    if (!ex) {
28        for (i = uxi << 12; i >> 63 == 0; ex--, i <<= 1)
29            ;
30        uxi <<= -ex + 1;
31    } else {
32        uxi &= -1ULL >> 12;
33        uxi |= 1ULL << 52;
34    }
35    if (!ey) {
36        for (i = uy.i << 12; i >> 63 == 0; ey--, i <<= 1)
37            ;
38        uy.i <<= -ey + 1;
39    } else {
40        uy.i &= -1ULL >> 12;
41        uy.i |= 1ULL << 52;
42    }
43
44    /* x mod y */
45    for (; ex > ey; ex--) {
46        i = uxi - uy.i;
47        if (i >> 63 == 0) {
48            if (i == 0)
49                return 0 * x;
50            uxi = i;
51        }
52        uxi <<= 1;
53    }
54    i = uxi - uy.i;
55    if (i >> 63 == 0) {
56        if (i == 0)
57            return 0 * x;
58        uxi = i;
59    }
60    for (; uxi >> 52 == 0; uxi <<= 1, ex--)
61        ;
62
63    /* scale result */
64    if (ex > 0) {
65        uxi -= 1ULL << 52;
66        uxi |= (uint64_t)ex << 52;
67    } else {
68        uxi >>= -ex + 1;
69    }
70    uxi |= (uint64_t)sx << 63;
71    ux.i = uxi;
72    return ux.f;
73}
74