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