libm-test.inc: Computing ulps near FP_ZERO.
Carlos O'Donell
carlos@redhat.com
Tue Apr 9 01:04:00 GMT 2013
On 04/07/2013 08:28 PM, Brooks Moses wrote:
> Carlos O'Donell wrote, at 4/7/2013 3:58 PM:
>> Excellent discussion.
>
> Thanks! Glad it was helpful.
>
>> At this point I'm in agreement with you.
>>
>> The likely course of action is to do the following:
>>
>> (a) Adjust tests that rely on accurate input to generate accurate
>> output, such as the cos(pi/2) test. These tests will have their
>> expected result adjusted by the error involved in converting
>> their input argument.
>>
>> (b) Enhance algorithms to detect zero branches. In particular fix
>> cpow to detect zero branches. We have a microbenchmark framework
>> now so I could run the new algorithms through that to see if
>> the function is faster or slower.
>>
>> Comments?
>
> FWIW, that sounds reasonable to me, but I'm just the peanut gallery. :)
>
> I'd be interested in knowing how the microbenchmarks come out. My guess
> is that adding check for zero-branches (which you'd implement just by
> checking for imag(z)==0, and then do other stuff only if that's true)
> will be negligible cost on most architectures with branch prediction.
> And my suggested alternate implementation might be faster, because
> you're avoiding some trigonometric calls.
I'll post results as further patches.
>> What about the cases where we can't easily fix the algorithm to
>> detect the zero branches? What do we use for ulp(0x0.0p0) then? :-)
>
> I figure those fall into three categories:
>
> * You're getting zero through something that amounts to subtracting
> two "equal" quantities (in the above cases, that was a pi/2 input
> and a pi/2 used to condition the input, basically), and the
> expected error is based on ulp(input), not ulp(output).
>
> * A zero input propagates through to a zero output. Here, you ought
> to get exactly zero in the output. This is what happened in the
> algorithm I wrote out for w^z; it doesn't actually detect zeros in
> a branching sense, but will preserve them nonetheless.
>
> * A zero output results from underflow. Here you also ought to get
> exactly zero.
>
> Are there any cases in the testsuite that are exceptions to this
> categorization?
I don't think so.
> I think anything where you can't get an exact zero is going to end up
> being like the w^z-with-the-current-implementation case -- where the
> error is inherently on the order of ulp(1) or so, and the best the test
> can do is document why -- rather than giving you an error on the order
> of the smallest floating-point numbers.
The following is the only testsuite failures I see after fixing up
the definition of ulp to match what is expected from a mathematical
perspective.
It reveals quite clearly, that cos, sincos, and cpow have some
terrible errors.
Regenerating ULPs for /home/carlos/build/glibc/math/test-idouble
testing double (inline functions)
Failure: Test: cos (pi/2) == PI_ERROR/2
Result:
is: 6.12323399573676603587e-17 0x1.1a62633145c070000000p-54
should be: 1.11022302462515654042e-16 0x1.00000000000000000000p-53
difference: 4.97899625051479936837e-17 0x1.cb3b399d747f20000000p-55
ulp(x) : 2.46519032881566189191e-32 0x1.00000000000000000000p-105
ulp : 2019720827359740.5000
max.ulp : 0.0000
Maximal error of `cos'
is : 2019720827359741 ulp
accepted: 2 ulp
Failure: Test: sincos (pi/2, &sin_res, &cos_res) puts PI_ERROR/2 in cos_res
Result:
is: 6.12323399573676603587e-17 0x1.1a62633145c070000000p-54
should be: 1.11022302462515654042e-16 0x1.00000000000000000000p-53
difference: 4.97899625051479936837e-17 0x1.cb3b399d747f20000000p-55
ulp(x) : 2.46519032881566189191e-32 0x1.00000000000000000000p-105
ulp : 2019720827359740.5000
max.ulp : 0.0000
Maximal error of `sincos'
is : 2019720827359741 ulp
accepted: 1 ulp
Failure: Test: Imaginary part of: cpow (e + 0 i, 0 + 2 * M_PIl i) == 1.0 - PI_ERROR i
Result:
is: -2.44929359829470641435e-16 -0x1.1a62633145c070000000p-52
should be: -2.22044604925031308085e-16 -0x1.00000000000000000000p-52
difference: 2.28847549044393333500e-17 0x1.a62633145c0700000000p-56
ulp(x) : 4.93038065763132378382e-32 0x1.00000000000000000000p-104
ulp : 464157972651015.0000
max.ulp : 0.0000
Maximal error of real part of: cpow
is : 2 ulp
accepted: 2 ulp
Maximal error of imaginary part of: cpow
is : 464157972651015 ulp
accepted: 2 ulp
Test suite completed:
7881 test cases plus 7253 tests for exception flags executed.
6 errors occurred.
Cheers,
Carlos.
-------------- next part --------------
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 78d2107..816682a 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,53 @@ 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))
+
+#ifdef TEST_FLOAT
+# define PI_ERROR (- ulp (M_PIl)/2)
+#endif
+#if defined TEST_DOUBLE || defined TEST_LDOUBLE
+# define PI_ERROR (ulp (M_PIl)/2)
+#endif
+
+/* Returns the size of an ulp for VALUE. */
+static FLOAT
+ulp (FLOAT value)
+{
+ FLOAT ulp;
+
+ switch (fpclassify (value))
+ {
+ case FP_ZERO:
+ /* If we compute the distance to the next FP it will be the same as the
+ value of the smallest subnormal number. This would yield a tiny value
+ the would mean that any answer not near zero would have a huge number
+ of ULPs. Instead we decide that the next FP is that of the nearest
+ normal value i.e. not subnormal. 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. The next normal value is going to
+ be 2^(1 - MAX_EXP). This lines up with the Java implementation of ulp. */
+ ulp = FUNC(ldexp) (1.0, 1 - MAX_EXP);
+ break;
+
+ case FP_NORMAL:
+ ulp = FUNC(ldexp) (1.0, FUNC(ilogb) (value) - MANT_DIG);
+ break;
+
+ case FP_SUBNORMAL:
+ /* The next closest subnormal value is a constant distance away. */
+ ulp = FUNC(ldexp) (1.0, 1 - (MAX_EXP + 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 +594,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 +622,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 +652,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);
}
}
@@ -5330,7 +5363,7 @@ cos_test (void)
TEST_f_f (cos, M_PI_6l * 2.0, 0.5);
TEST_f_f (cos, M_PI_6l * 4.0, -0.5);
- TEST_f_f (cos, M_PI_2l, 0);
+ TEST_f_f (cos, M_PI_2l, PI_ERROR/2);
TEST_f_f (cos, 0.75L, 0.731688868873820886311838753000084544L);
@@ -5651,7 +5684,7 @@ 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);
- TEST_cc_c (cpow, M_El, 0, 0, 2 * M_PIl, 1.0, 0.0);
+ TEST_cc_c (cpow, M_El, 0, 0, 2 * M_PIl, 1.0, -PI_ERROR);
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);
@@ -12133,8 +12166,8 @@ sincos_test (void)
TEST_extra (sincos, plus_infty, qnan_value, qnan_value, INVALID_EXCEPTION);
TEST_extra (sincos, minus_infty, qnan_value, qnan_value, INVALID_EXCEPTION);
TEST_extra (sincos, qnan_value, qnan_value, qnan_value);
+ TEST_extra (sincos, M_PI_2l, 1, PI_ERROR/2);
- TEST_extra (sincos, M_PI_2l, 1, 0);
TEST_extra (sincos, M_PI_6l, 0.5, 0.86602540378443864676372317075293616L);
TEST_extra (sincos, M_PI_6l*2.0, 0.86602540378443864676372317075293616L, 0.5);
TEST_extra (sincos, 0.75L, 0.681638760023334166733241952779893935L, 0.731688868873820886311838753000084544L);
@@ -13045,31 +13078,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 subnormal value... */
+ ulps = ulp (0x0.0p0);
+ if (fpclassify (ulps) != FP_NORMAL)
+ {
+ fprintf (stderr, "ulp (0x0.0p0) is not FP_NORMAL!\n");
+ exit (EXIT_FAILURE);
+ }
+ /* ... and that it is the smallest possible subnormal. */
+ value = FUNC (ldexp) (1.0, 1 - MAX_EXP) ;
+ if (ulps != value)
+ {
+ fprintf (stderr, "ulp (0x0.0p0) is not 2^(1 - MAX_EXP)\n");
+ exit (EXIT_FAILURE);
+ }
}
-#endif
int
main (int argc, char **argv)
@@ -13120,9 +13149,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: */
More information about the Libc-alpha
mailing list