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