Correct IBM long double nextafterl.

Fix for values near a power of two, and some tidies.

	[BZ #16739]
	* sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c (__nextafterl): Correct
	output when value is near a power of two.  Use int64_t for lx and
	remove casts.  Use decimal rather than hex exponent constants.
	Don't use long double multiplication when double will suffice.
	* math/libm-test.inc (nextafter_test_data): Add tests.
	* NEWS: Add 16739 and 16786 to bug list.

diff --git a/ChangeLog b/ChangeLog
index 7c737abaf8cb1124111efc7e1da149511b896770..f623e8b79a30142077112a66bc6844dc2424c8e4 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,13 @@
+2014-04-02  Alan Modra  <amodra@gmail.com>
+
+       [BZ #16739]
+       * sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c (__nextafterl): Correct
+       output when value is near a power of two.  Use int64_t for lx and
+       remove casts.  Use decimal rather than hex exponent constants.
+       Don't use long double multiplication when double will suffice.
+       * math/libm-test.inc (nextafter_test_data): Add tests.
+       * NEWS: Add 16739 and 16786 to bug list.
+
 2014-04-02  Alan Modra  <amodra@gmail.com>
 
        * sysdeps/powerpc/powerpc64/power7/memrchr.S: Correct stream hint.

diff --git a/NEWS b/NEWS
index 31723768ebc24fe66e8227e3d3eda4d032b599b5..d959fb33779799bc3271b6617257e47900621091 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -13,7 +13,7 @@ Version 2.20
   16357, 16362, 16447, 16532, 16545, 16574, 16599, 16600, 16609, 16610,
   16611, 16613, 16623, 16632, 16634, 16639, 16642, 16648, 16649, 16670,
   16674, 16677, 16680, 16683, 16689, 16695, 16701, 16706, 16707, 16712,
-  16713, 16714, 16731, 16743, 16758, 16759, 16760, 16770.
+  16713, 16714, 16731, 16739, 16743, 16758, 16759, 16760, 16770, 16786.
 
 * Running the testsuite no longer terminates as soon as a test fails.
   Instead, a file tests.sum (xtests.sum from "make xcheck") is generated,

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 19194f63e9e07577f2501783ba5cebf5ed360c4c..967b679e9407ab5f0b88ff046fc92387933c1eaf 100644 (file)
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -8302,6 +8302,14 @@ static const struct test_ff_f_data nextafter_test_data[] =
     // XXX Enable once gcc is fixed.
     //TEST_ff_f (nextafter, 0x0.00000040000000000000p-16385L, -0.1L, 0x0.0000003ffffffff00000p-16385L),
 #endif
+#if defined TEST_LDOUBLE && LDBL_MANT_DIG == 106
+    TEST_ff_f (nextafter, 1.0L, -10.0L, 1.0L-0x1p-106L, NO_EXCEPTION),
+    TEST_ff_f (nextafter, 1.0L, 10.0L, 1.0L+0x1p-105L, NO_EXCEPTION),
+    TEST_ff_f (nextafter, 1.0L-0x1p-106L, 10.0L, 1.0L, NO_EXCEPTION),
+    TEST_ff_f (nextafter, -1.0L, -10.0L, -1.0L-0x1p-105L, NO_EXCEPTION),
+    TEST_ff_f (nextafter, -1.0L, 10.0L, -1.0L+0x1p-106L, NO_EXCEPTION),
+    TEST_ff_f (nextafter, -1.0L+0x1p-106L, -10.0L, -1.0L, NO_EXCEPTION),
+#endif
 
     /* XXX We need the hexadecimal FP number representation here for further
        tests.  */

diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c b/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c
index 30b1540a8849a4c9d9b07b852b22950fb4faec49..bf57cb89d62e92af65cd08dbabef33c42bf5faff 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c
@@ -30,8 +30,7 @@ static char rcsid[] = "$NetBSD: $";
 
 long double __nextafterl(long double x, long double y)
 {
-       int64_t hx,hy,ihx,ihy;
-       uint64_t lx;
+       int64_t hx, hy, ihx, ihy, lx;
        double xhi, xlo, yhi;
 
        ldbl_unpack (x, &xhi, &xlo);
@@ -79,19 +78,28 @@ long double __nextafterl(long double x, long double y)
              u = math_opt_barrier (x);
              x -= __LDBL_DENORM_MIN__;
              if (ihx < 0x0360000000000000LL
-                 || (hx > 0 && (int64_t) lx <= 0)
-                 || (hx < 0 && (int64_t) lx > 1)) {
+                 || (hx > 0 && lx <= 0)
+                 || (hx < 0 && lx > 1)) {
                u = u * u;
                math_force_eval (u);            /* raise underflow flag */
              }
              return x;
            }
-           if (ihx < 0x06a0000000000000LL) { /* ulp will denormal */
-             INSERT_WORDS64 (yhi, hx & (0x7ffLL<<52));
-             u = yhi;
-             u *= 0x1.0000000000000p-105L;
+           /* If the high double is an exact power of two and the low
+              double is the opposite sign, then 1ulp is one less than
+              what we might determine from the high double.  Similarly
+              if X is an exact power of two, and positive, because
+              making it a little smaller will result in the exponent
+              decreasing by one and normalisation of the mantissa.   */
+           if ((hx & 0x000fffffffffffffLL) == 0
+               && ((lx != 0 && (hx ^ lx) < 0)
+                   || (lx == 0 && hx >= 0)))
+             ihx -= 1LL << 52;
+           if (ihx < (106LL << 52)) { /* ulp will denormal */
+             INSERT_WORDS64 (yhi, ihx & (0x7ffLL<<52));
+             u = yhi * 0x1p-105;
            } else {
-             INSERT_WORDS64 (yhi, (hx & (0x7ffLL<<52))-(0x069LL<<52));
+             INSERT_WORDS64 (yhi, (ihx & (0x7ffLL<<52))-(105LL<<52));
              u = yhi;
            }
            return x - u;
@@ -109,8 +117,8 @@ long double __nextafterl(long double x, long double y)
              u = math_opt_barrier (x);
              x += __LDBL_DENORM_MIN__;
              if (ihx < 0x0360000000000000LL
-                 || (hx > 0 && (int64_t) lx < 0 && lx != 0x8000000000000001LL)
-                 || (hx < 0 && (int64_t) lx >= 0)) {
+                 || (hx > 0 && lx < 0 && lx != 0x8000000000000001LL)
+                 || (hx < 0 && lx >= 0)) {
                u = u * u;
                math_force_eval (u);            /* raise underflow flag */
              }
@@ -118,12 +126,21 @@ long double __nextafterl(long double x, long double y)
                x = -0.0L;
              return x;
            }
-           if (ihx < 0x06a0000000000000LL) { /* ulp will denormal */
-             INSERT_WORDS64 (yhi, hx & (0x7ffLL<<52));
-             u = yhi;
-             u *= 0x1.0000000000000p-105L;
+           /* If the high double is an exact power of two and the low
+              double is the opposite sign, then 1ulp is one less than
+              what we might determine from the high double.  Similarly
+              if X is an exact power of two, and negative, because
+              making it a little larger will result in the exponent
+              decreasing by one and normalisation of the mantissa.   */
+           if ((hx & 0x000fffffffffffffLL) == 0
+               && ((lx != 0 && (hx ^ lx) < 0)
+                   || (lx == 0 && hx < 0)))
+             ihx -= 1LL << 52;
+           if (ihx < (106LL << 52)) { /* ulp will denormal */
+             INSERT_WORDS64 (yhi, ihx & (0x7ffLL<<52));
+             u = yhi * 0x1p-105;
            } else {
-             INSERT_WORDS64 (yhi, (hx & (0x7ffLL<<52))-(0x069LL<<52));
+             INSERT_WORDS64 (yhi, (ihx & (0x7ffLL<<52))-(105LL<<52));
              u = yhi;
            }
            return x + u;