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