This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[PATCH] libm-test.inc: Correctly implement ulp().
- From: "Carlos O'Donell" <carlos at redhat dot com>
- To: "Joseph S. Myers" <joseph at codesourcery dot com>, Andreas Jaeger <aj at suse dot com>, GNU C Library <libc-alpha at sourceware dot org>
- Date: Thu, 16 May 2013 18:10:02 -0400
- Subject: [PATCH] libm-test.inc: Correctly implement ulp().
Joseph, Andreas,
The following patch fixes the master TODO entry:
~~~
When libm-test.inc computes ulps for a case where the
expected result was 0, it computes them as if the expected
result was 1, but should compute them more strictly, as
if the expected result was subnormal.
~~~
The current value used for ulp around zero is just wrong,
and this test fixes it such that ulp(0) is the smallest subnormal
value nearest to zero, which makes the most sense for testing
values near zero.
The patch causes 1 regressions for each of the types on x86-64,
all having to do with " cpow (e + 0 i, 0 + 2 * M_PIl i) == 1.0 + 0.0 i".
e.g.
~~~
Regenerating ULPs for /home/carlos/build/glibc/math/test-double
testing double (without inline functions)
Failure: Test: Imaginary part of: cpow (e + 0 i, 0 + 2 * M_PIl i) == 1.0 + 0.0 i
Result:
is: -2.44929359829470641435e-16 -0x1.1a62633145c070000000p-52
should be: 0.00000000000000000000e+00 0x0.00000000000000000000p+0
difference: 2.44929359829470641435e-16 0x1.1a62633145c070000000p-52
ulp(x) : 4.94065645841246544177e-324 0x0.00000000000010000000p-1022
ulp : 49574254330601946645624577457449570499924092730262502276699108050758815961872456534680590923878394894568106529394830647426366788862820389065085884559356492887832038112641869699009951764897429803306592466752526966543116836508033509283588356790340998671355914085615325047212736424084369897623125485610193649664.0000
max.ulp : 2.0000
Maximal error of real part of: cpow
is : 2 ulp
accepted: 2 ulp
Maximal error of imaginary part of: cpow
is : 49574254330601946645624577457449570499924092730262502276699108050758815961872456534680590923878394894568106529394830647426366788862820389065085884559356492887832038112641869699009951764897429803306592466752526966543116836508033509283588356790340998671355914085615325047212736424084369897623125485610193649664 ulp
accepted: 2 ulp
Test suite completed:
8044 test cases plus 7416 tests for exception flags executed.
2 errors occurred.
~~~
I thought that perhaps I could use the same strategy as I had
used with sincos to fix this failure, but I can't. Even correcting
the answer for the approximation of e and 2 * M_PIl still results
in an imaginary answer that is not near enough to zero.
Instead as Joseph has suggested in the past I've disabled this
test and indicated that it is because of BZ #14473:
Bug 14473 - Inaccurate CPOW* function on all machines.
http://sourceware.org/bugzilla/show_bug.cgi?id=14473
No regressions on x86_64.
OK to checkin?
Cheers,
Carlos.
v2
- Disable failing cpow test.
- Enhance check_ulp to use nextafter for 1, 2, and 100 ulp.
2013-05-16 Carlos O'Donell <carlos@redhat.com>
* math/libm-test.inc (MAX_EXP): Define.
(ULPDIFF): Define.
(ulp): New function.
(check_float_internal): Use ULPDIFF, and
print ulp(expected) in output.
(cpow_test): Disable failing test.
(check_ulp): Test ulp() implemetnation.
(main): Call check_ulp before starting tests.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 3a7acb8..9887d46 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -270,6 +270,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))
static void
init_max_error (void)
@@ -589,6 +591,44 @@ test_errno (const char *test_name, int errno_value, int exceptions)
test_single_errno (test_name, errno_value, ERANGE, "ERANGE");
}
+/* Returns the number of ulps that GIVEN is away from EXPECTED. */
+#define ULPDIFF(given, expected) \
+ (FUNC(fabs) ((given) - (expected)) / ulp (expected))
+
+/* Returns the size of an ulp for VALUE. */
+static FLOAT
+ulp (FLOAT value)
+{
+ FLOAT ulp;
+
+ switch (fpclassify (value))
+ {
+ case FP_ZERO:
+ /* We compute the distance to the next FP which is the same as the
+ value of the smallest subnormal number. Previously we used
+ 2^(-MANT_DIG) which is too large a value to be useful. Note that we
+ can't use ilogb(0), since that isn't a valid thing to do. As a point
+ of comparison Java's ulp returns the next normal value e.g.
+ 2^(1 - MAX_EXP) for ulp(0), but that is not what we want for
+ glibc. */
+ /* 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));
+ break;
+
+ case FP_NORMAL:
+ ulp = FUNC(ldexp) (1.0, FUNC(ilogb) (value) - MANT_DIG);
+ break;
+
+ default:
+ /* It should never happen. */
+ abort ();
+ break;
+ }
+ return ulp;
+}
+
static void
check_float_internal (const char *test_name, FLOAT computed, FLOAT expected,
FLOAT max_ulp, int exceptions,
@@ -597,7 +637,7 @@ check_float_internal (const char *test_name, FLOAT computed, FLOAT expected,
int ok = 0;
int print_diff = 0;
FLOAT diff = 0;
- FLOAT ulp = 0;
+ FLOAT ulps = 0;
int errno_value = errno;
test_exceptions (test_name, exceptions);
@@ -627,37 +667,19 @@ check_float_internal (const char *test_name, FLOAT computed, FLOAT expected,
else
{
diff = FUNC(fabs) (computed - expected);
- switch (fpclassify (expected))
- {
- case FP_ZERO:
- /* ilogb (0) isn't allowed. */
- ulp = diff / FUNC(ldexp) (1.0, - MANT_DIG);
- break;
- case FP_NORMAL:
- ulp = diff / FUNC(ldexp) (1.0, FUNC(ilogb) (expected) - MANT_DIG);
- break;
- case FP_SUBNORMAL:
- /* 1ulp for a subnormal value, shifted by MANT_DIG, is the
- least normal value. */
- ulp = (FUNC(ldexp) (diff, MANT_DIG) / min_value);
- break;
- default:
- /* It should never happen. */
- abort ();
- break;
- }
- set_max_error (ulp, curr_max_error);
+ ulps = ULPDIFF (computed, expected);
+ set_max_error (ulps, curr_max_error);
print_diff = 1;
if ((exceptions & IGNORE_ZERO_INF_SIGN) == 0
&& computed == 0.0 && expected == 0.0
&& signbit(computed) != signbit (expected))
ok = 0;
- else if (ulp <= 0.5 || (ulp <= max_ulp && !ignore_max_ulp))
+ else if (ulps <= 0.5 || (ulps <= max_ulp && !ignore_max_ulp))
ok = 1;
else
{
ok = 0;
- print_ulps (test_name, ulp);
+ print_ulps (test_name, ulps);
}
}
@@ -675,7 +697,9 @@ check_float_internal (const char *test_name, FLOAT computed, FLOAT expected,
{
printf (" difference: % .20" PRINTF_EXPR " % .20" PRINTF_XEXPR
"\n", diff, diff);
- printf (" ulp : % .4" PRINTF_NEXPR "\n", ulp);
+ printf (" ulp(x) : % .20" PRINTF_EXPR " % .20" PRINTF_XEXPR
+ "\n", ulp (expected), ulp (expected));
+ printf (" ulp : % .4" PRINTF_NEXPR "\n", ulps);
printf (" max.ulp : % .4" PRINTF_NEXPR "\n", max_ulp);
}
}
@@ -6991,7 +7015,11 @@ cpow_test (void)
TEST_cc_c (cpow, 1, 0, 0, 0, 1.0, 0.0);
TEST_cc_c (cpow, 2, 0, 10, 0, 1024.0, 0.0);
+ /* Disabled until we fix BZ #14473 */
+#if 0
TEST_cc_c (cpow, M_El, 0, 0, 2 * M_PIl, 1.0, 0.0);
+#endif
+
TEST_cc_c (cpow, 2, 3, 4, 0, -119.0, -120.0);
TEST_cc_c (cpow, qnan_value, qnan_value, qnan_value, qnan_value, qnan_value, qnan_value);
@@ -14910,31 +14938,54 @@ parse_opt (int key, char *arg, struct argp_state *state)
return 0;
}
-#if 0
-/* function to check our ulp calculation. */
+/* Verify that our ulp () implementation is behaving as expected
+ or abort. */
void
check_ulp (void)
{
- int i;
-
- FLOAT u, diff, ulp;
- /* This gives one ulp. */
- u = FUNC(nextafter) (10, 20);
- check_equal (10.0, u, 1, &diff, &ulp);
- printf ("One ulp: % .4" PRINTF_NEXPR "\n", ulp);
-
- /* This gives one more ulp. */
- u = FUNC(nextafter) (u, 20);
- check_equal (10.0, u, 2, &diff, &ulp);
- printf ("two ulp: % .4" PRINTF_NEXPR "\n", ulp);
-
- /* And now calculate 100 ulp. */
- for (i = 2; i < 100; i++)
- u = FUNC(nextafter) (u, 20);
- check_equal (10.0, u, 100, &diff, &ulp);
- printf ("100 ulp: % .4" PRINTF_NEXPR "\n", ulp);
+ FLOAT ulps, value;
+ int i;
+ /* Check ulp of zero is a subnormal value... */
+ ulps = ulp (0x0.0p0);
+ if (fpclassify (ulps) != FP_SUBNORMAL)
+ {
+ fprintf (stderr, "ulp (0x0.0p0) is not FP_SUBNORMAL!\n");
+ exit (EXIT_FAILURE);
+ }
+ /* Check that the ulp of one is a normal value... */
+ ulps = ulp (1.0L);
+ if (fpclassify (ulps) != FP_NORMAL)
+ {
+ fprintf (stderr, "ulp (1.0L) is not FP_NORMAL\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. */
+ value = FUNC(nextafter) (10, 20);
+ ulps = ULPDIFF (value, 10);
+ if (ulps > 1.0L)
+ {
+ fprintf (stderr, "The value of ulp (10+1ulp,10) is greater than 1 ulp.\n");
+ exit (EXIT_FAILURE);
+ }
+ /* This gives one more ulp. */
+ value = FUNC(nextafter) (value, 20);
+ ulps = ULPDIFF (value, 10);
+ if (ulps > 2.0L)
+ {
+ fprintf (stderr, "The value of ulp (10+2ulp,10) is greater than 2 ulp.\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)
+ {
+ fprintf (stderr, "The value of ulp (10+100ulp,10) is greater than 100 ulp.\n");
+ exit (EXIT_FAILURE);
+ }
}
-#endif
int
main (int argc, char **argv)
@@ -14985,9 +15036,7 @@ main (int argc, char **argv)
initialize ();
printf (TEST_MSG);
-#if 0
check_ulp ();
-#endif
/* Keep the tests a wee bit ordered (according to ISO C99). */
/* Classification macros: */
---
Cheers,
Carlos.