Re: [PATCH] BZ#15536: Fix ulp(0x0.0p0) for 128-bit IBM long double.

["Carlos O'Donell" <carlos at redhat dot com>]

On 05/31/2013 03:07 PM, Joseph S. Myers wrote:
> On Fri, 31 May 2013, Carlos O'Donell wrote:
> 
>> @@ -269,8 +269,8 @@ static FLOAT max_error, real_max_error, imag_max_error;
>>  
>>  #define MANT_DIG CHOOSE ((LDBL_MANT_DIG-1), (DBL_MANT_DIG-1), (FLT_MANT_DIG-1),  \
>>                          (LDBL_MANT_DIG-1), (DBL_MANT_DIG-1), (FLT_MANT_DIG-1))
>> -#define MAX_EXP CHOOSE ((LDBL_MAX_EXP-1), (DBL_MAX_EXP-1), (FLT_MAX_EXP-1),    \
>> -                       (LDBL_MAX_EXP-1), (DBL_MAX_EXP-1), (FLT_MAX_EXP-1))
>> +#define MIN_EXP CHOOSE ((LDBL_MIN_EXP), (DBL_MIN_EXP), (FLT_MIN_EXP),  \
>> +                       (LDBL_MIN_EXP), (DBL_MIN_EXP), (FLT_MIN_EXP))
>>  
>>  /* Compare KEY (a string, with the name of a test or a function) with
>>     ULP (a pointer to a struct ulp_data structure), returning a value
>> @@ -680,7 +680,7 @@ ulp (FLOAT value)
>>         /* Fall through...  */
>>        case FP_SUBNORMAL:
>>          /* The next closest subnormal value is a constant distance away.  */
>> -       ulp = FUNC(ldexp) (1.0, 1 - (MAX_EXP + MANT_DIG));
>> +       ulp = FUNC(ldexp) (1.0, MIN_EXP - MANT_DIG);
> 
> You have an off-by-one error here, so this would give e.g. 0x1p-1073 not 
> 0x1p-1074 for double.
> 
> (The confusing semantics of the ISO C <float.h> macros don't particularly 
> help, but it's worse when the macros in libm-test.inc mysteriously 
> subtract 1 so their semantics are different from the ISO C semantics, and 
> even worse when some macros but not others subtract 1.)

Yes, you're right, I forgot MANT_DIG subtracts one. I really wanted
*_MANT_DIG unmodified for this particular solution.

It is annoying that the semantics are different, I'll make a note that
this should really be cleaned up to use the semantics as defined for
ISO C, not some arbitrary semantics (or if we do, name them something
different like MIN_EXP_M1).

I've updated check_ulp to ensure this doesn't happen again. The first
implementation of check ulp didn't have any strict checking for ulp
smaller than expected, just ulp larger than expected. I've fixed it
to check for value +/- 1ulp.

I verified that with the check_ulp changes the test fails to run
because ulp is wrong. After fixing MIN_EXP to fix the off-by-one 
error the check_ulp sanity test passes and the test runs without
error again (modulo imprecise results from complex functions
for 128-bit IBM long double which is BZ #15396).

Passes on x86-64 and Power64.

OK to checkin?

v2
- Use DBL_MANT_DIG.
v3
- Use MIN_EXP - MANT_DIG.
v4
- Fix check_ulp to detect error.
- Define MIN_EXP as `*_MIN_EXP - 1`.

2013-05-31  Carlos O'Donell  <carlos@redhat.com>

	[BZ #15536]
	* math/libm-test.inc (MAX_EXP): Remove
	(MIN_EXP): Define.
	(ulp): Use MIN_EXP - MANT_DIG.
	(check_ulp): Verify subnormal ulps. Only allow expected range of +/-1ulp.

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 7a6bf09..6870d96 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -269,8 +269,8 @@ static FLOAT max_error, real_max_error, imag_max_error;
 
 #define MANT_DIG CHOOSE ((LDBL_MANT_DIG-1), (DBL_MANT_DIG-1), (FLT_MANT_DIG-1),  \
                         (LDBL_MANT_DIG-1), (DBL_MANT_DIG-1), (FLT_MANT_DIG-1))
-#define MAX_EXP CHOOSE ((LDBL_MAX_EXP-1), (DBL_MAX_EXP-1), (FLT_MAX_EXP-1),    \
-                       (LDBL_MAX_EXP-1), (DBL_MAX_EXP-1), (FLT_MAX_EXP-1))
+#define MIN_EXP CHOOSE ((LDBL_MIN_EXP-1), (DBL_MIN_EXP-1), (FLT_MIN_EXP-1),    \
+                       (LDBL_MIN_EXP-1), (DBL_MIN_EXP-1), (FLT_MIN_EXP-1))
 
 /* Compare KEY (a string, with the name of a test or a function) with
    ULP (a pointer to a struct ulp_data structure), returning a value
@@ -680,7 +680,7 @@ ulp (FLOAT value)
        /* Fall through...  */
       case FP_SUBNORMAL:
         /* The next closest subnormal value is a constant distance away.  */
-       ulp = FUNC(ldexp) (1.0, 1 - (MAX_EXP + MANT_DIG));
+       ulp = FUNC(ldexp) (1.0, MIN_EXP - MANT_DIG);
        break;
 
       case FP_NORMAL:
@@ -14583,7 +14583,7 @@ parse_opt (int key, char *arg, struct argp_state *state)
 void
 check_ulp (void)
 {
-   FLOAT ulps, value;
+   FLOAT ulps, ulpx, value;
    int i;
    /* Check ulp of zero is a subnormal value...  */
    ulps = ulp (0x0.0p0);
@@ -14599,30 +14599,45 @@ check_ulp (void)
        fprintf (stderr, "ulp (1.0L) is not FP_NORMAL\n");
        exit (EXIT_FAILURE);
      }
+
+   /* Compute the next subnormal value using nextafter to validate ulp.
+      We allow +/- 1 ulp around the represented value.  */
+   value = FUNC(nextafter) (0, 1);
+   ulps = ULPDIFF (value, 0);
+   ulpx = ulp (1.0L);
+   if (ulps < (1.0L - ulpx) || ulps > (1.0L + ulpx))
+     {
+       fprintf (stderr, "Value outside of 1 +/- 1ulp.\n");
+       exit (EXIT_FAILURE);
+     }
    /* Compute the nearest representable number from 10 towards 20.
-      The result is 10 + 1ulp.  We use this to check the ulp function.  */
+      The result is 10 + 1ulp.  We use this to check the ulp function.
+      We allow +/- 1 ulp around the represented value.  */
    value = FUNC(nextafter) (10, 20);
    ulps = ULPDIFF (value, 10);
-   if (ulps > 1.0L)
+   ulpx = ulp (1.0L);
+   if (ulps < (1.0L - ulpx) || ulps > (1.0L + ulpx))
      {
-       fprintf (stderr, "The value of ulp (10+1ulp,10) is greater than 1 ulp.\n");
+       fprintf (stderr, "Value outside of 1 +/- 1ulp.\n");
        exit (EXIT_FAILURE);
      }
    /* This gives one more ulp.  */
    value = FUNC(nextafter) (value, 20);
    ulps = ULPDIFF (value, 10);
-   if (ulps > 2.0L)
+   ulpx = ulp (2.0L);
+   if (ulps < (2.0L - ulpx) || ulps > (2.0L + ulpx))
      {
-       fprintf (stderr, "The value of ulp (10+2ulp,10) is greater than 2 ulp.\n");
+       fprintf (stderr, "Value outside of 2 +/- 1ulp.\n");
        exit (EXIT_FAILURE);
      }
    /* And now calculate 100 ulp.  */
    for (i = 2; i < 100; i++)
      value = FUNC(nextafter) (value, 20);
    ulps = ULPDIFF (value, 10);
-   if (ulps > 100.0L)
+   ulpx = ulp (100.0L);
+   if (ulps < (100.0L - ulpx) || ulps > (100.0L + ulpx))
      {
-       fprintf (stderr, "The value of ulp (10+100ulp,10) is greater than 100 ulp.\n");
+       fprintf (stderr, "Value outside of 100 +/- 1ulp.\n");
        exit (EXIT_FAILURE);
      }
 }
---

Cheers,
Carlos.