OLD | NEW |
1 #include "libm.h" | 1 #include "libm.h" |
2 | 2 |
3 #if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1 | 3 #if FLT_EVAL_METHOD == 0 || FLT_EVAL_METHOD == 1 |
4 #define EPS DBL_EPSILON | 4 #define EPS DBL_EPSILON |
5 #elif FLT_EVAL_METHOD==2 | 5 #elif FLT_EVAL_METHOD == 2 |
6 #define EPS LDBL_EPSILON | 6 #define EPS LDBL_EPSILON |
7 #endif | 7 #endif |
8 static const double_t toint = 1/EPS; | 8 static const double_t toint = 1 / EPS; |
9 | 9 |
10 double floor(double x) | 10 double floor(double x) { |
11 { | 11 union { |
12 » union {double f; uint64_t i;} u = {x}; | 12 double f; |
13 » int e = u.i >> 52 & 0x7ff; | 13 uint64_t i; |
14 » double_t y; | 14 } u = {x}; |
| 15 int e = u.i >> 52 & 0x7ff; |
| 16 double_t y; |
15 | 17 |
16 » if (e >= 0x3ff+52 || x == 0) | 18 if (e >= 0x3ff + 52 || x == 0) |
17 » » return x; | 19 return x; |
18 » /* y = int(x) - x, where int(x) is an integer neighbor of x */ | 20 /* y = int(x) - x, where int(x) is an integer neighbor of x */ |
19 » if (u.i >> 63) | 21 if (u.i >> 63) |
20 » » y = x - toint + toint - x; | 22 y = x - toint + toint - x; |
21 » else | 23 else |
22 » » y = x + toint - toint - x; | 24 y = x + toint - toint - x; |
23 » /* special case because of non-nearest rounding modes */ | 25 /* special case because of non-nearest rounding modes */ |
24 » if (e <= 0x3ff-1) { | 26 if (e <= 0x3ff - 1) { |
25 » » FORCE_EVAL(y); | 27 FORCE_EVAL(y); |
26 » » return u.i >> 63 ? -1 : 0; | 28 return u.i >> 63 ? -1 : 0; |
27 » } | 29 } |
28 » if (y > 0) | 30 if (y > 0) |
29 » » return x + y - 1; | 31 return x + y - 1; |
30 » return x + y; | 32 return x + y; |
31 } | 33 } |
OLD | NEW |