This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[WIP] Fixing ulp near zero.
- From: "Carlos O'Donell" <carlos at redhat dot com>
- To: "Joseph S. Myers" <joseph at codesourcery dot com>, GNU C Library <libc-alpha at sourceware dot org>
- Date: Thu, 25 Apr 2013 00:07:57 -0400
- Subject: [WIP] Fixing ulp near zero.
Joseph,
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 patch causes 2 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.
~~~
Where the expected answer of 0.0i is naive and is a similar
problem to rounding pi/2.
I plan on fixing the cpow result per your recommendations
from the cos and sincos changes such that the ulp error is
at least smaller.
Comments?
Cheers,
Carlos.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index eb9fa71..2127f88 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -248,6 +248,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)
@@ -536,6 +538,43 @@ test_exceptions (const char *test_name, int exception)
feclearexcept (FE_ALL_EXCEPT);
}
+/* 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,
@@ -545,7 +584,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;
test_exceptions (test_name, exceptions);
if (issignaling (computed) && issignaling (expected))
@@ -573,37 +612,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);
}
}
@@ -621,7 +642,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);
}
}
@@ -13169,31 +13192,27 @@ 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;
+ /* Check ulp of zero is a normal value... */
+ ulps = ulp (0x0.0p0);
+ if (fpclassify (ulps) != FP_SUBNORMAL)
+ {
+ fprintf (stderr, "ulp (0x0.0p0) is not FP_SUBNORMAL!\n");
+ exit (EXIT_FAILURE);
+ }
+ /* ... and that it is the smallest possible subnormal. */
+ value = FUNC (ldexp) (1.0, 1 - MAX_EXP - MANT_DIG) ;
+ if (ulps != value)
+ {
+ fprintf (stderr, "ulp (0x0.0p0) is not 2^(1 - MAX_EXP - MANT_DIG)\n");
+ exit (EXIT_FAILURE);
+ }
}
-#endif
int
main (int argc, char **argv)
@@ -13244,9 +13263,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: */
---