diff mbox

[PowerPC64] Correct IBM long double nextafterl

Message ID 20140325022618.GF18201@bubble.grove.modra.org
State New
Headers show

Commit Message

Alan Modra March 25, 2014, 2:26 a.m. UTC
Fix for values near a power of two, and some tidies.

The added tests do cause some new fails, because IBM long double isn't
supposed to support rounding modes other than round-to-nearest.
I believe powerpc should have a math-tests.h like arm and mips that
limits the tested rounding modes.

	[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.

Comments

Joseph Myers March 25, 2014, 12:13 p.m. UTC | #1
On Tue, 25 Mar 2014, Alan Modra wrote:

> +    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),

You mean NO_INEXACT_EXCEPTION (NO_EXCEPTION makes no assertion about 
whether or not "inexact" is present; "inexact" has a different default 
from the other exceptions, following ISO C Annex F).
Alan Modra March 25, 2014, 1:15 p.m. UTC | #2
On Tue, Mar 25, 2014 at 12:13:08PM +0000, Joseph S. Myers wrote:
> On Tue, 25 Mar 2014, Alan Modra wrote:
> 
> > +    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),
> 
> You mean NO_INEXACT_EXCEPTION (NO_EXCEPTION makes no assertion about 
> whether or not "inexact" is present; "inexact" has a different default 
> from the other exceptions, following ISO C Annex F).

I deliberately did not choose NO_INEXACT_EXCEPTION, because nextafterl
on IBM long double *does* set FE_INEXACT on some of these tests.  The
reason is the "x + u" or "x - u" operations inside nextafterl which
use __gcc_qadd.  As you know, __gcc_qadd spuriously sets FE_INEXACT.
Also, as the description of IBM long double in the PowerPC64 ABI
states, the "Extended precision" format

 * Does not support the IEEE status flags for overflow, underflow, and 
   other conditions.  These flag have no meaning in this format.
Adhemerval Zanella March 27, 2014, 8:11 p.m. UTC | #3
On 25-03-2014 10:15, Alan Modra wrote:
> On Tue, Mar 25, 2014 at 12:13:08PM +0000, Joseph S. Myers wrote:
>> On Tue, 25 Mar 2014, Alan Modra wrote:
>>
>>> +    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),
>> You mean NO_INEXACT_EXCEPTION (NO_EXCEPTION makes no assertion about 
>> whether or not "inexact" is present; "inexact" has a different default 
>> from the other exceptions, following ISO C Annex F).
> I deliberately did not choose NO_INEXACT_EXCEPTION, because nextafterl
> on IBM long double *does* set FE_INEXACT on some of these tests.  The
> reason is the "x + u" or "x - u" operations inside nextafterl which
> use __gcc_qadd.  As you know, __gcc_qadd spuriously sets FE_INEXACT.
> Also, as the description of IBM long double in the PowerPC64 ABI
> states, the "Extended precision" format
>
>  * Does not support the IEEE status flags for overflow, underflow, and 
>    other conditions.  These flag have no meaning in this format.
>
If Joseph any more concerns, I'm ok with this patch.
Joseph Myers March 27, 2014, 9:55 p.m. UTC | #4
On Thu, 27 Mar 2014, Adhemerval Zanella wrote:

> If Joseph any more concerns, I'm ok with this patch.

I have no further concerns with this patch.
diff mbox

Patch

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 5e50f0e..e97b18a 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -8312,6 +8312,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 30b1540..bf57cb8 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);
@@ -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;