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