From e266b71770050a4d0cb276f4afea1c5b05215184 Mon Sep 17 00:00:00 2001 From: Alan Modra Date: Wed, 2 Apr 2014 13:46:19 +1030 Subject: [PATCH] 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. --- ChangeLog | 10 +++++ NEWS | 2 +- math/libm-test.inc | 8 ++++ sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c | 49 +++++++++++++++------- 4 files changed, 52 insertions(+), 17 deletions(-) diff --git a/ChangeLog b/ChangeLog index fe00e0c3159..5bf3e2667ab 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,13 @@ +2014-04-02 Alan Modra + + [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 * sysdeps/powerpc/powerpc64/power7/memrchr.S: Correct stream hint. diff --git a/NEWS b/NEWS index a647e5d736a..71f6fb6618d 100644 --- 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. Version 2.19 diff --git a/math/libm-test.inc b/math/libm-test.inc index 1c9b97a8e77..cf07069d8c3 100644 --- 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 c050944c0c2..7b09927e26e 100644 --- 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; -- 2.47.2