Fix nextafter/nexttoward overflow in directed rounding modes

src/s_nextafter.c

@@ -64,7 +64,7 @@ nextafter(double x, double y)
 	    }
 	}
 	hy = hx&0x7ff00000;
-	if(hy>=0x7ff00000) return x+x;	/* overflow  */
+	if(hy>=0x7ff00000) t = x+x;	/* raise overflow flag */
 	if(hy<0x00100000) {		/* underflow */
 	    t = x*x;
 	    if(t!=x) {		/* raise underflow flag */

src/s_nextafterf.c

@@ -54,7 +54,7 @@ nextafterf(float x, float y)
 	    }
 	}
 	hy = hx&0x7f800000;
-	if(hy>=0x7f800000) return x+x;	/* overflow  */
+	if(hy>=0x7f800000) t = x+x;	/* raise overflow flag */
 	if(hy<0x00800000) {		/* underflow */
 	    t = x*x;
 	    if(t!=x) {		/* raise underflow flag */

src/s_nextafterl.c

@@ -67,7 +67,7 @@ nextafterl(long double x, long double y)
 		    ux.bits.exp += 1;
 	    }
 	}
-	if(ux.bits.exp==0x7fff) return x+x;	/* overflow  */
+	if(ux.bits.exp==0x7fff) t = x+x;	/* raise overflow flag */
 	if(ux.bits.exp==0) {			/* underflow */
 	    mask_nbit_l(ux);
 	    t = ux.e * ux.e;

src/s_nexttoward.c

@@ -59,7 +59,7 @@ nexttoward(double x, long double y)
 	    if(lx==0) hx += 1;
 	}
 	ix = hx&0x7ff00000;
-	if(ix>=0x7ff00000) return x+x;	/* overflow  */
+	if(ix>=0x7ff00000) t = x+x;	/* raise overflow flag */
 	if(ix<0x00100000) {		/* underflow */
 	    t = x*x;
 	    if(t!=x) {		/* raise underflow flag */

src/s_nexttowardf.c

@@ -47,7 +47,7 @@ nexttowardf(float x, long double y)
 	else					/* x += ulp */
 	    hx += 1;
 	ix = hx&0x7f800000;
-	if(ix>=0x7f800000) return x+x;	/* overflow  */
+	if(ix>=0x7f800000) t = x+x;	/* raise overflow flag */
 	if(ix<0x00800000) {		/* underflow */
 	    t = x*x;
 	    if(t!=x) {		/* raise underflow flag */

test/regression/test-361.c

@@ -0,0 +1,59 @@
+/*
+ * Regression test for issue #361: nextafter/nexttoward must return signed
+ * infinity at the finite endpoints in every rounding mode, while raising
+ * overflow and inexact (C Annex F.10.8.3-4).
+ */
+#include <float.h>
+#include <openlibm_fenv.h>
+#include <openlibm_math.h>
+
+#include "regress-util.h"
+
+#define CHECK_ENDPOINT(expr) do { \
+	long double result; \
+	int flags; \
+	CHECK(feclearexcept(FE_ALL_EXCEPT) == 0); \
+	result = (expr); \
+	flags = fetestexcept(FE_ALL_EXCEPT); \
+	CHECK(isinf(result)); \
+	CHECK(!!signbit(result) == (sign < 0)); \
+	CHECK((flags & (FE_OVERFLOW | FE_INEXACT)) == \
+	    (FE_OVERFLOW | FE_INEXACT)); \
+	CHECK(fegetround() == modes[m]); \
+} while (0)
+
+int
+main(void)
+{
+	const int modes[] = {
+		FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO
+	};
+	int saved_round = fegetround();
+	unsigned int m;
+	int sign;
+
+	for (m = 0; m < sizeof(modes) / sizeof(modes[0]); m++) {
+		CHECK(fesetround(modes[m]) == 0);
+		for (sign = -1; sign <= 1; sign += 2) {
+			CHECK_ENDPOINT(nextafter(sign * DBL_MAX,
+			    sign * INFINITY));
+			CHECK_ENDPOINT(nextafterf(sign * FLT_MAX,
+			    sign * INFINITY));
+#ifdef OLM_LONG_DOUBLE
+			CHECK_ENDPOINT(nexttowardf(sign * FLT_MAX,
+			    sign * (long double)INFINITY));
+/* Extended long double routines are built only on x86 and aarch64. */
+#if LDBL_MANT_DIG == 53 || LDBL_MANT_DIG == 64 || defined(__aarch64__)
+			CHECK_ENDPOINT(nextafterl(sign * LDBL_MAX,
+			    sign * (long double)INFINITY));
+			CHECK_ENDPOINT(nexttoward(sign * DBL_MAX,
+			    sign * (long double)INFINITY));
+			CHECK_ENDPOINT(nexttowardl(sign * LDBL_MAX,
+			    sign * (long double)INFINITY));
+#endif
+#endif
+		}
+	}
+	CHECK(fesetround(saved_round) == 0);
+	return 0;
+}