Fix dbl-64 exp overflow/underflow in non-default rounding modes (bug 16284)
Siddhesh Poyarekar
siddhesh@redhat.com
Mon Mar 24 09:15:00 GMT 2014
On Fri, Mar 21, 2014 at 11:24:26PM +0000, Joseph S. Myers wrote:
> The dbl-64 version of exp needs round-to-nearest mode for its internal
> computations, but that has the consequence of inappropriate
> overflowing and underflowing results in other rounding modes. This
> patch fixes this by recomputing the relevant results in cases where
> the round-to-nearest result overflows to infinity or underflows to
> zero (most of the diffs are actually just consequent reindentation).
> Tests are enabled in all rounding modes for complex functions using
> exp - but not for cexp because it turns out there are bugs causing
> spurious underflows for cexp for some tests, which will need to be
> fixed separately (I suspect ccos ccosh csin csinh ctan ctanh have
> similar bugs, just not shown by the present set of test inputs).
>
> Tested x86_64 and x86 and ulps updated accordingly.
>
> (auto-libm-test-out diffs omitted below.)
>
> 2014-03-21 Joseph Myers <joseph@codesourcery.com>
>
> [BZ #16284]
> * sysdeps/ieee754/dbl-64/e_exp.c (__ieee754_exp): Use original
> rounding mode to recompute results that overflow to infinity or
> underflow to zero.
> * math/auto-libm-test-in: Don't mark tests as expected to fail for
> bug 16284.
> * math/auto-libm-test-out: Regenerated.
> * math/libm-test.inc (ccos_test): Use ALL_RM_TEST.
> (ccosh_test): Likewise.
> (csin_test_data): Use plus_oflow.
> (csin_test): Use ALL_RM_TEST.
> (csinh_test_data): Use plus_oflow.
> (csinh_test): Use ALL_RM_TEST.
> * sysdeps/i386/fpu/libm-test-ulps: Update.
> * sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
>
Looks good to me.
> diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
> index fafe96f..5975631 100644
> --- a/math/auto-libm-test-in
> +++ b/math/auto-libm-test-in
> @@ -833,16 +833,14 @@ exp 0.75
> exp 50.0
> exp 88.72269439697265625
> exp 709.75
> -# Bug 16284: results on directed rounding may be incorrect.
> # GCC bug 59666: results on directed rounding may be incorrect.
> -exp 1000.0 xfail-rounding:dbl-64 xfail-rounding:ldbl-128ibm
> -exp 710 xfail-rounding:dbl-64 xfail-rounding:ldbl-128ibm
> +exp 1000.0 xfail-rounding:ldbl-128ibm
> +exp 710 xfail-rounding:ldbl-128ibm
> exp -1234
> -# Bug 16284: results on directed rounding may be incorrect.
> # GCC bug 59666: results on directed rounding may be incorrect.
> -exp 0x2.c679d1f73f0fb628p+8 xfail-rounding:dbl-64 xfail-rounding:ldbl-128ibm
> -exp 1e5 xfail-rounding:dbl-64 xfail-rounding:ldbl-128ibm
> -exp max xfail-rounding:dbl-64 xfail-rounding:ldbl-128ibm
> +exp 0x2.c679d1f73f0fb628p+8 xfail-rounding:ldbl-128ibm
> +exp 1e5 xfail-rounding:ldbl-128ibm
> +exp max xfail-rounding:ldbl-128ibm
> exp -7.4444006192138124e+02
> exp -0x1.75f113c30b1c8p+9
> exp -max
> @@ -856,27 +854,22 @@ exp10 36
> exp10 -36
> exp10 305
> exp10 -305
> -# Bug 16284: results on directed rounding may be incorrect.
> # GCC bug 59666: results on directed rounding may be incorrect.
> -exp10 4932 xfail-rounding:flt-32 xfail-rounding:ldbl-128ibm
> +exp10 4932 xfail-rounding:ldbl-128ibm
> # Bug 16361: underflow exception may be misssing
> exp10 -4932 missing-underflow:ldbl-96-intel:x86 missing-underflow:ldbl-96-intel:x86_64
> -# Bug 16284: results on directed rounding may be incorrect.
> # GCC bug 59666: results on directed rounding may be incorrect.
> -exp10 1e5 xfail-rounding:flt-32 xfail-rounding:ldbl-128ibm
> +exp10 1e5 xfail-rounding:ldbl-128ibm
> exp10 -1e5
> -# Bug 16284: results on directed rounding may be incorrect.
> # GCC bug 59666: results on directed rounding may be incorrect.
> -exp10 1e6 xfail-rounding:flt-32 xfail-rounding:ldbl-128ibm
> +exp10 1e6 xfail-rounding:ldbl-128ibm
> exp10 -1e6
> -# Bug 16284: results on directed rounding may be incorrect.
> # GCC bug 59666: results on directed rounding may be incorrect.
> -exp10 max xfail-rounding:flt-32 xfail-rounding:ldbl-128ibm
> +exp10 max xfail-rounding:ldbl-128ibm
> exp10 -max
> exp10 0.75
> -# Bug 16284: results on directed rounding may be incorrect.
> # GCC bug 59666: results on directed rounding may be incorrect.
> -exp10 0x1.348e45573a1dd72cp+8 xfail-rounding:flt-32 xfail-rounding:dbl-64 xfail-rounding:ldbl-128ibm
> +exp10 0x1.348e45573a1dd72cp+8 xfail-rounding:ldbl-128ibm
>
> exp2 0
> exp2 -0
> diff --git a/math/libm-test.inc b/math/libm-test.inc
> index 5e50f0e..a8ebecd 100644
> --- a/math/libm-test.inc
> +++ b/math/libm-test.inc
> @@ -5820,9 +5820,7 @@ static const struct test_c_c_data ccos_test_data[] =
> static void
> ccos_test (void)
> {
> - START (ccos, 0);
> - RUN_TEST_LOOP_c_c (ccos, ccos_test_data, );
> - END_COMPLEX;
> + ALL_RM_TEST (ccos, 0, ccos_test_data, RUN_TEST_LOOP_c_c, END_COMPLEX);
> }
>
>
> @@ -5879,9 +5877,7 @@ static const struct test_c_c_data ccosh_test_data[] =
> static void
> ccosh_test (void)
> {
> - START (ccosh, 0);
> - RUN_TEST_LOOP_c_c (ccosh, ccosh_test_data, );
> - END_COMPLEX;
> + ALL_RM_TEST (ccosh, 0, ccosh_test_data, RUN_TEST_LOOP_c_c, END_COMPLEX);
> }
>
>
> @@ -6467,15 +6463,15 @@ static const struct test_c_c_data csin_test_data[] =
> #endif
>
> #ifdef TEST_FLOAT
> - TEST_c_c (csin, 0x1p-149, 180, 1.043535896672617552965983803453927655332e33L, plus_infty, OVERFLOW_EXCEPTION),
> + TEST_c_c (csin, 0x1p-149, 180, 1.043535896672617552965983803453927655332e33L, plus_oflow, OVERFLOW_EXCEPTION),
> #endif
>
> #if defined TEST_DOUBLE || (defined TEST_LDOUBLE && LDBL_MAX_EXP == 1024)
> - TEST_c_c (csin, 0x1p-1074, 1440, 5.981479269486130556466515778180916082415e301L, plus_infty, OVERFLOW_EXCEPTION),
> + TEST_c_c (csin, 0x1p-1074, 1440, 5.981479269486130556466515778180916082415e301L, plus_oflow, OVERFLOW_EXCEPTION),
> #endif
>
> #if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
> - TEST_c_c (csin, 0x1p-16434L, 22730, 1.217853148905605987081057582351152052687e4924L, plus_infty, OVERFLOW_EXCEPTION),
> + TEST_c_c (csin, 0x1p-16434L, 22730, 1.217853148905605987081057582351152052687e4924L, plus_oflow, OVERFLOW_EXCEPTION),
> #endif
>
> TEST_c_c (csin, min_subnorm_value, min_value, min_subnorm_value, min_value, UNDERFLOW_EXCEPTION),
> @@ -6485,9 +6481,7 @@ static const struct test_c_c_data csin_test_data[] =
> static void
> csin_test (void)
> {
> - START (csin, 0);
> - RUN_TEST_LOOP_c_c (csin, csin_test_data, );
> - END_COMPLEX;
> + ALL_RM_TEST (csin, 0, csin_test_data, RUN_TEST_LOOP_c_c, END_COMPLEX);
> }
>
>
> @@ -6566,15 +6560,15 @@ static const struct test_c_c_data csinh_test_data[] =
> #endif
>
> #ifdef TEST_FLOAT
> - TEST_c_c (csinh, 180, 0x1p-149, plus_infty, 1.043535896672617552965983803453927655332e33L, OVERFLOW_EXCEPTION),
> + TEST_c_c (csinh, 180, 0x1p-149, plus_oflow, 1.043535896672617552965983803453927655332e33L, OVERFLOW_EXCEPTION),
> #endif
>
> #if defined TEST_DOUBLE || (defined TEST_LDOUBLE && LDBL_MAX_EXP == 1024)
> - TEST_c_c (csinh, 1440, 0x1p-1074, plus_infty, 5.981479269486130556466515778180916082415e301L, OVERFLOW_EXCEPTION),
> + TEST_c_c (csinh, 1440, 0x1p-1074, plus_oflow, 5.981479269486130556466515778180916082415e301L, OVERFLOW_EXCEPTION),
> #endif
>
> #if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
> - TEST_c_c (csinh, 22730, 0x1p-16434L, plus_infty, 1.217853148905605987081057582351152052687e4924L, OVERFLOW_EXCEPTION),
> + TEST_c_c (csinh, 22730, 0x1p-16434L, plus_oflow, 1.217853148905605987081057582351152052687e4924L, OVERFLOW_EXCEPTION),
> #endif
>
> TEST_c_c (csinh, min_subnorm_value, min_value, min_subnorm_value, min_value, UNDERFLOW_EXCEPTION),
> @@ -6584,9 +6578,7 @@ static const struct test_c_c_data csinh_test_data[] =
> static void
> csinh_test (void)
> {
> - START (csinh, 0);
> - RUN_TEST_LOOP_c_c (csinh, csinh_test_data, );
> - END_COMPLEX;
> + ALL_RM_TEST (csinh, 0, csinh_test_data, RUN_TEST_LOOP_c_c, END_COMPLEX);
> }
>
>
> diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
> index b7b2e12..1885be7 100644
> --- a/sysdeps/i386/fpu/libm-test-ulps
> +++ b/sysdeps/i386/fpu/libm-test-ulps
> @@ -437,6 +437,54 @@ ifloat: 1
> ildouble: 1
> ldouble: 1
>
> +Function: Real part of "ccos_downward":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "ccos_downward":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "ccos_towardzero":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "ccos_towardzero":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "ccos_upward":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 2
> +ldouble: 2
> +
> +Function: Imaginary part of "ccos_upward":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 2
> +ldouble: 2
> +
> Function: Real part of "ccosh":
> double: 1
> float: 1
> @@ -451,6 +499,54 @@ ifloat: 1
> ildouble: 1
> ldouble: 1
>
> +Function: Real part of "ccosh_downward":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "ccosh_downward":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "ccosh_towardzero":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "ccosh_towardzero":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "ccosh_upward":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 2
> +ldouble: 2
> +
> +Function: Imaginary part of "ccosh_upward":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 2
> +ldouble: 2
> +
> Function: Real part of "cexp":
> double: 1
> float: 1
> @@ -582,6 +678,54 @@ float: 1
> idouble: 1
> ifloat: 1
>
> +Function: Real part of "csin_downward":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csin_downward":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "csin_towardzero":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csin_towardzero":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "csin_upward":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csin_upward":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> Function: Real part of "csinh":
> double: 1
> float: 1
> @@ -596,6 +740,54 @@ float: 1
> idouble: 1
> ifloat: 1
>
> +Function: Real part of "csinh_downward":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csinh_downward":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "csinh_towardzero":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csinh_towardzero":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "csinh_upward":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csinh_upward":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> Function: Real part of "csqrt":
> double: 1
> idouble: 1
> diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c
> index 100fc39..cd060ce 100644
> --- a/sysdeps/ieee754/dbl-64/e_exp.c
> +++ b/sysdeps/ieee754/dbl-64/e_exp.c
> @@ -58,168 +58,178 @@ __ieee754_exp (double x)
> int4 i, j, m, n, ex;
> double retval;
>
> - SET_RESTORE_ROUND (FE_TONEAREST);
> + {
> + SET_RESTORE_ROUND (FE_TONEAREST);
>
> - junk1.x = x;
> - m = junk1.i[HIGH_HALF];
> - n = m & hugeint;
> + junk1.x = x;
> + m = junk1.i[HIGH_HALF];
> + n = m & hugeint;
>
> - if (n > smallint && n < bigint)
> - {
> - y = x * log2e.x + three51.x;
> - bexp = y - three51.x; /* multiply the result by 2**bexp */
> + if (n > smallint && n < bigint)
> + {
> + y = x * log2e.x + three51.x;
> + bexp = y - three51.x; /* multiply the result by 2**bexp */
>
> - junk1.x = y;
> + junk1.x = y;
>
> - eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */
> - t = x - bexp * ln_two1.x;
> + eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */
> + t = x - bexp * ln_two1.x;
>
> - y = t + three33.x;
> - base = y - three33.x; /* t rounded to a multiple of 2**-18 */
> - junk2.x = y;
> - del = (t - base) - eps; /* x = bexp*ln(2) + base + del */
> - eps = del + del * del * (p3.x * del + p2.x);
> + y = t + three33.x;
> + base = y - three33.x; /* t rounded to a multiple of 2**-18 */
> + junk2.x = y;
> + del = (t - base) - eps; /* x = bexp*ln(2) + base + del */
> + eps = del + del * del * (p3.x * del + p2.x);
>
> - binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
> + binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
>
> - i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
> - j = (junk2.i[LOW_HALF] & 511) << 1;
> + i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
> + j = (junk2.i[LOW_HALF] & 511) << 1;
>
> - al = coar.x[i] * fine.x[j];
> - bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
> - + coar.x[i + 1] * fine.x[j + 1]);
> + al = coar.x[i] * fine.x[j];
> + bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
> + + coar.x[i + 1] * fine.x[j + 1]);
>
> - rem = (bet + bet * eps) + al * eps;
> - res = al + rem;
> - cor = (al - res) + rem;
> - if (res == (res + cor * err_0))
> - {
> - retval = res * binexp.x;
> - goto ret;
> - }
> - else
> - {
> - retval = __slowexp (x);
> - goto ret;
> - } /*if error is over bound */
> - }
> + rem = (bet + bet * eps) + al * eps;
> + res = al + rem;
> + cor = (al - res) + rem;
> + if (res == (res + cor * err_0))
> + {
> + retval = res * binexp.x;
> + goto ret;
> + }
> + else
> + {
> + retval = __slowexp (x);
> + goto ret;
> + } /*if error is over bound */
> + }
>
> - if (n <= smallint)
> - {
> - retval = 1.0;
> - goto ret;
> - }
> + if (n <= smallint)
> + {
> + retval = 1.0;
> + goto ret;
> + }
>
> - if (n >= badint)
> - {
> - if (n > infint)
> - {
> - retval = x + x;
> - goto ret;
> - } /* x is NaN */
> - if (n < infint)
> - {
> - retval = (x > 0) ? (hhuge * hhuge) : (tiny * tiny);
> - goto ret;
> - }
> - /* x is finite, cause either overflow or underflow */
> - if (junk1.i[LOW_HALF] != 0)
> - {
> - retval = x + x;
> - goto ret;
> - } /* x is NaN */
> - retval = (x > 0) ? inf.x : zero; /* |x| = inf; return either inf or 0 */
> - goto ret;
> - }
> + if (n >= badint)
> + {
> + if (n > infint)
> + {
> + retval = x + x;
> + goto ret;
> + } /* x is NaN */
> + if (n < infint)
> + {
> + if (x > 0)
> + goto ret_huge;
> + else
> + goto ret_tiny;
OK.
> + }
> + /* x is finite, cause either overflow or underflow */
> + if (junk1.i[LOW_HALF] != 0)
> + {
> + retval = x + x;
> + goto ret;
> + } /* x is NaN */
> + retval = (x > 0) ? inf.x : zero; /* |x| = inf; return either inf or 0 */
> + goto ret;
> + }
>
> - y = x * log2e.x + three51.x;
> - bexp = y - three51.x;
> - junk1.x = y;
> - eps = bexp * ln_two2.x;
> - t = x - bexp * ln_two1.x;
> - y = t + three33.x;
> - base = y - three33.x;
> - junk2.x = y;
> - del = (t - base) - eps;
> - eps = del + del * del * (p3.x * del + p2.x);
> - i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
> - j = (junk2.i[LOW_HALF] & 511) << 1;
> - al = coar.x[i] * fine.x[j];
> - bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
> - + coar.x[i + 1] * fine.x[j + 1]);
> - rem = (bet + bet * eps) + al * eps;
> - res = al + rem;
> - cor = (al - res) + rem;
> - if (m >> 31)
> - {
> - ex = junk1.i[LOW_HALF];
> - if (res < 1.0)
> - {
> - res += res;
> - cor += cor;
> - ex -= 1;
> - }
> - if (ex >= -1022)
> - {
> - binexp.i[HIGH_HALF] = (1023 + ex) << 20;
> - if (res == (res + cor * err_0))
> - {
> - retval = res * binexp.x;
> - goto ret;
> - }
> - else
> - {
> - retval = __slowexp (x);
> - goto check_uflow_ret;
> - } /*if error is over bound */
> - }
> - ex = -(1022 + ex);
> - binexp.i[HIGH_HALF] = (1023 - ex) << 20;
> - res *= binexp.x;
> - cor *= binexp.x;
> - eps = 1.0000000001 + err_0 * binexp.x;
> - t = 1.0 + res;
> - y = ((1.0 - t) + res) + cor;
> - res = t + y;
> - cor = (t - res) + y;
> - if (res == (res + eps * cor))
> - {
> - binexp.i[HIGH_HALF] = 0x00100000;
> - retval = (res - 1.0) * binexp.x;
> - goto check_uflow_ret;
> - }
> - else
> - {
> - retval = __slowexp (x);
> - goto check_uflow_ret;
> - } /* if error is over bound */
> - check_uflow_ret:
> - if (retval < DBL_MIN)
> - {
> + y = x * log2e.x + three51.x;
> + bexp = y - three51.x;
> + junk1.x = y;
> + eps = bexp * ln_two2.x;
> + t = x - bexp * ln_two1.x;
> + y = t + three33.x;
> + base = y - three33.x;
> + junk2.x = y;
> + del = (t - base) - eps;
> + eps = del + del * del * (p3.x * del + p2.x);
> + i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
> + j = (junk2.i[LOW_HALF] & 511) << 1;
> + al = coar.x[i] * fine.x[j];
> + bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
> + + coar.x[i + 1] * fine.x[j + 1]);
> + rem = (bet + bet * eps) + al * eps;
> + res = al + rem;
> + cor = (al - res) + rem;
> + if (m >> 31)
> + {
> + ex = junk1.i[LOW_HALF];
> + if (res < 1.0)
> + {
> + res += res;
> + cor += cor;
> + ex -= 1;
> + }
> + if (ex >= -1022)
> + {
> + binexp.i[HIGH_HALF] = (1023 + ex) << 20;
> + if (res == (res + cor * err_0))
> + {
> + retval = res * binexp.x;
> + goto ret;
> + }
> + else
> + {
> + retval = __slowexp (x);
> + goto check_uflow_ret;
> + } /*if error is over bound */
> + }
> + ex = -(1022 + ex);
> + binexp.i[HIGH_HALF] = (1023 - ex) << 20;
> + res *= binexp.x;
> + cor *= binexp.x;
> + eps = 1.0000000001 + err_0 * binexp.x;
> + t = 1.0 + res;
> + y = ((1.0 - t) + res) + cor;
> + res = t + y;
> + cor = (t - res) + y;
> + if (res == (res + eps * cor))
> + {
> + binexp.i[HIGH_HALF] = 0x00100000;
> + retval = (res - 1.0) * binexp.x;
> + goto check_uflow_ret;
> + }
> + else
> + {
> + retval = __slowexp (x);
> + goto check_uflow_ret;
This is probably a separate cleanup, but the two jumps above could be
removed.
> + } /* if error is over bound */
> + check_uflow_ret:
> + if (retval < DBL_MIN)
> + {
> #if FLT_EVAL_METHOD != 0
> - volatile
> + volatile
> #endif
> - double force_underflow = tiny * tiny;
> - math_force_eval (force_underflow);
> - }
> - goto ret;
> - }
> - else
> - {
> - binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
> - if (res == (res + cor * err_0))
> - {
> + double force_underflow = tiny * tiny;
> + math_force_eval (force_underflow);
> + }
> + if (retval == 0)
> + goto ret_tiny;
OK.
> + goto ret;
> + }
> + else
> + {
> + binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
> + if (res == (res + cor * err_0))
> retval = res * binexp.x * t256.x;
> - goto ret;
> - }
> - else
> - {
> + else
> retval = __slowexp (x);
> + if (__isinf (retval))
> + goto ret_huge;
OK.
> + else
> goto ret;
> - }
> - }
> + }
> + }
> ret:
> return retval;
> +
> + ret_huge:
> + return hhuge * hhuge;
> +
> + ret_tiny:
> + return tiny * tiny;
OK.
> }
> #ifndef __ieee754_exp
> strong_alias (__ieee754_exp, __exp_finite)
> diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
> index 301eaa6..670f2da 100644
> --- a/sysdeps/x86_64/fpu/libm-test-ulps
> +++ b/sysdeps/x86_64/fpu/libm-test-ulps
> @@ -468,6 +468,54 @@ ifloat: 1
> ildouble: 1
> ldouble: 1
>
> +Function: Real part of "ccos_downward":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "ccos_downward":
> +double: 2
> +float: 3
> +idouble: 2
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "ccos_towardzero":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "ccos_towardzero":
> +double: 2
> +float: 3
> +idouble: 2
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "ccos_upward":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 2
> +ldouble: 2
> +
> +Function: Imaginary part of "ccos_upward":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 2
> +ldouble: 2
> +
> Function: Real part of "ccosh":
> double: 1
> float: 1
> @@ -482,6 +530,54 @@ ifloat: 1
> ildouble: 1
> ldouble: 1
>
> +Function: Real part of "ccosh_downward":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "ccosh_downward":
> +double: 2
> +float: 3
> +idouble: 2
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "ccosh_towardzero":
> +double: 1
> +float: 3
> +idouble: 1
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "ccosh_towardzero":
> +double: 2
> +float: 3
> +idouble: 2
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "ccosh_upward":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 2
> +ldouble: 2
> +
> +Function: Imaginary part of "ccosh_upward":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 2
> +ldouble: 2
> +
> Function: Real part of "cexp":
> double: 2
> float: 1
> @@ -616,6 +712,54 @@ ifloat: 1
> ildouble: 1
> ldouble: 1
>
> +Function: Real part of "csin_downward":
> +double: 2
> +float: 3
> +idouble: 2
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csin_downward":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "csin_towardzero":
> +double: 2
> +float: 3
> +idouble: 2
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csin_towardzero":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "csin_upward":
> +double: 1
> +float: 3
> +idouble: 1
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csin_upward":
> +double: 1
> +float: 3
> +idouble: 1
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> Function: Real part of "csinh":
> float: 1
> ifloat: 1
> @@ -628,6 +772,54 @@ float: 1
> idouble: 1
> ifloat: 1
>
> +Function: Real part of "csinh_downward":
> +double: 1
> +float: 2
> +idouble: 1
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csinh_downward":
> +double: 2
> +float: 3
> +idouble: 2
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "csinh_towardzero":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csinh_towardzero":
> +double: 2
> +float: 3
> +idouble: 2
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Real part of "csinh_upward":
> +double: 1
> +float: 3
> +idouble: 1
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> +Function: Imaginary part of "csinh_upward":
> +double: 2
> +float: 3
> +idouble: 2
> +ifloat: 3
> +ildouble: 3
> +ldouble: 3
> +
> Function: Real part of "csqrt":
> double: 1
> float: 1
>
> --
> Joseph S. Myers
> joseph@codesourcery.com
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 473 bytes
Desc: not available
URL: <http://sourceware.org/pipermail/libc-alpha/attachments/20140324/1cd6fa04/attachment.sig>
More information about the Libc-alpha
mailing list