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.

Backport of b0abbc21034f0e5edc49023d8fda0616173faf17.

diff --git a/ChangeLog b/ChangeLog
index fe00e0c31592cd88a857fb75c35ae8c95816a3d9..5bf3e2667ab39aaa0c1ce89ebd2abf1b004f4afe 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 a647e5d736a1cb88466f1eefef021b62cf2b67b9..71f6fb6618ddf2c9cea468e8ae222d8858a6a207 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -9,7 +9,7 @@ Version 2.19.1
 
 * The following bugs are resolved with this release:
 
-  16545, 16683, 16689, 16701, 16706, 16707.
+  16545, 16683, 16689, 16701, 16706, 16707, 16739.
 \f
 Version 2.19
 

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 1c9b97a8e77907c2468800c39b6360176aca51c5..cf07069d8c39d162475b7912fe9fa942cad9d574 100644 (file)
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -10547,6 +10547,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 c050944c0c270134fcd5b73d8f0eb8709fbad249..7b09927e26ef7448ca87f4b6f36cb91a51b4fd82 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);
@@ -76,19 +75,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;
@@ -103,8 +111,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 */
              }
@@ -112,12 +120,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;