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