1 #include <math.h>
2 #include <stdint.h>
3 
fmod(double x,double y)4 double 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