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: */