OLD | NEW |
1 #include "libm.h" | 1 #include "libm.h" |
2 | 2 |
3 double nextafter(double x, double y) | 3 double nextafter(double x, double y) { |
4 { | 4 union { |
5 » union {double f; uint64_t i;} ux={x}, uy={y}; | 5 double f; |
6 » uint64_t ax, ay; | 6 uint64_t i; |
7 » int e; | 7 } ux = {x}, uy = {y}; |
| 8 uint64_t ax, ay; |
| 9 int e; |
8 | 10 |
9 » if (isnan(x) || isnan(y)) | 11 if (isnan(x) || isnan(y)) |
10 » » return x + y; | 12 return x + y; |
11 » if (ux.i == uy.i) | 13 if (ux.i == uy.i) |
12 » » return y; | 14 return y; |
13 » ax = ux.i & -1ULL/2; | 15 ax = ux.i & -1ULL / 2; |
14 » ay = uy.i & -1ULL/2; | 16 ay = uy.i & -1ULL / 2; |
15 » if (ax == 0) { | 17 if (ax == 0) { |
16 » » if (ay == 0) | 18 if (ay == 0) |
17 » » » return y; | 19 return y; |
18 » » ux.i = (uy.i & 1ULL<<63) | 1; | 20 ux.i = (uy.i & 1ULL << 63) | 1; |
19 » } else if (ax > ay || ((ux.i ^ uy.i) & 1ULL<<63)) | 21 } else if (ax > ay || ((ux.i ^ uy.i) & 1ULL << 63)) |
20 » » ux.i--; | 22 ux.i--; |
21 » else | 23 else |
22 » » ux.i++; | 24 ux.i++; |
23 » e = ux.i >> 52 & 0x7ff; | 25 e = ux.i >> 52 & 0x7ff; |
24 » /* raise overflow if ux.f is infinite and x is finite */ | 26 /* raise overflow if ux.f is infinite and x is finite */ |
25 » if (e == 0x7ff) | 27 if (e == 0x7ff) |
26 » » FORCE_EVAL(x+x); | 28 FORCE_EVAL(x + x); |
27 » /* raise underflow if ux.f is subnormal or zero */ | 29 /* raise underflow if ux.f is subnormal or zero */ |
28 » if (e == 0) | 30 if (e == 0) |
29 » » FORCE_EVAL(x*x + ux.f*ux.f); | 31 FORCE_EVAL(x * x + ux.f * ux.f); |
30 » return ux.f; | 32 return ux.f; |
31 } | 33 } |
OLD | NEW |