ENH: add nextafterl implementation.

1 files changed, 54 insertions, 7 deletions

diff --git a/numpy/core/src/npymath/ieee754.c b/numpy/core/src/npymath/ieee754.c
index 902e8919b..c5b91c099 100644
--- a/numpy/core/src/npymath/ieee754.c
+++ b/numpy/core/src/npymath/ieee754.c
@@ -156,16 +156,63 @@ float npy_nextafterf(float x, float y)
 #endif
 
 #ifndef HAVE_NEXTAFTERL
-#if NPY_BITSOF_LONGDOUBLE == NPY_BITSOF_DOUBLE
 npy_longdouble npy_nextafterl(npy_longdouble x, npy_longdouble y)
 {
-    return (npy_longdouble)npy_nextafter((double)x, (double)y);
+    volatile npy_longdouble t;
+    union IEEEl2bits ux, uy;
+
+    ux.e = x;
+    uy.e = y;
+
+    if ((ux.bits.exp == 0x7fff &&
+	 ((ux.bits.manh & ~LDBL_NBIT) | ux.bits.manl) != 0) ||
+	(uy.bits.exp == 0x7fff &&
+	 ((uy.bits.manh & ~LDBL_NBIT) | uy.bits.manl) != 0)) {
+	return x + y;		/* x or y is nan */
+    }
+    if (x == y) {
+	return y;		/* x=y, return y */
+    }
+    if (x == 0.0) {
+	ux.bits.manh = 0;	/* return +-minsubnormal */
+	ux.bits.manl = 1;
+	ux.bits.sign = uy.bits.sign;
+	t = ux.e * ux.e;
+	if (t == ux.e) {
+	    return t;
+        } else {
+	    return ux.e;	/* raise underflow flag */
+        }
+    }
+    if (x > 0.0 ^ x < y) {	/* x -= ulp */
+	if (ux.bits.manl == 0) {
+	    if ((ux.bits.manh & ~LDBL_NBIT) == 0) {
+		ux.bits.exp -= 1;
+            }
+	    ux.bits.manh = (ux.bits.manh - 1) | (ux.bits.manh & LDBL_NBIT);
+	}
+	ux.bits.manl -= 1;
+    } else {			/* x += ulp */
+	ux.bits.manl += 1;
+	if (ux.bits.manl == 0) {
+	    ux.bits.manh = (ux.bits.manh + 1) | (ux.bits.manh & LDBL_NBIT);
+	    if ((ux.bits.manh & ~LDBL_NBIT) == 0) {
+		ux.bits.exp += 1;
+            }
+	}
+    }
+    if (ux.bits.exp == 0x7fff) {
+	return x + x;		/* overflow  */
+    }
+    if (ux.bits.exp == 0) {	/* underflow */
+	mask_nbit_l(ux);
+	t = ux.e * ux.e;
+	if (t != ux.e) {        /* raise underflow flag */
+	    return ux.e;
+        }
+    }
+    return ux.e;
 }
-#else
-/* long double is not standardized: we need to know the exact binary
- * representation for this platform */
-#error Needs nextafterl implementation for this platform
-#endif
 #endif
 
 /*