| #include <math.h> |
| #include <stdint.h> |
| |
| double fmod(double x, double y) { |
| union { |
| double f; |
| uint64_t i; |
| } ux = {x}, uy = {y}; |
| int ex = ux.i >> 52 & 0x7ff; |
| int ey = uy.i >> 52 & 0x7ff; |
| int sx = ux.i >> 63; |
| uint64_t i; |
| |
| /* in the followings uxi should be ux.i, but then gcc wrongly adds */ |
| /* float load/store to inner loops ruining performance and code size */ |
| uint64_t uxi = ux.i; |
| |
| if (uy.i << 1 == 0 || isnan(y) || ex == 0x7ff) |
| return (x * y) / (x * y); |
| if (uxi << 1 <= uy.i << 1) { |
| if (uxi << 1 == uy.i << 1) |
| return 0 * x; |
| return x; |
| } |
| |
| /* normalize x and y */ |
| if (!ex) { |
| for (i = uxi << 12; i >> 63 == 0; ex--, i <<= 1) |
| ; |
| uxi <<= -ex + 1; |
| } else { |
| uxi &= -1ULL >> 12; |
| uxi |= 1ULL << 52; |
| } |
| if (!ey) { |
| for (i = uy.i << 12; i >> 63 == 0; ey--, i <<= 1) |
| ; |
| uy.i <<= -ey + 1; |
| } else { |
| uy.i &= -1ULL >> 12; |
| uy.i |= 1ULL << 52; |
| } |
| |
| /* x mod y */ |
| for (; ex > ey; ex--) { |
| i = uxi - uy.i; |
| if (i >> 63 == 0) { |
| if (i == 0) |
| return 0 * x; |
| uxi = i; |
| } |
| uxi <<= 1; |
| } |
| i = uxi - uy.i; |
| if (i >> 63 == 0) { |
| if (i == 0) |
| return 0 * x; |
| uxi = i; |
| } |
| for (; uxi >> 52 == 0; uxi <<= 1, ex--) |
| ; |
| |
| /* scale result */ |
| if (ex > 0) { |
| uxi -= 1ULL << 52; |
| uxi |= (uint64_t)ex << 52; |
| } else { |
| uxi >>= -ex + 1; |
| } |
| uxi |= (uint64_t)sx << 63; |
| ux.i = uxi; |
| return ux.f; |
| } |