1 #include "libm.h"
2 
nextafter(double x,double y)3 double 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