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