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