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
43