| Message ID | Pine.LNX.4.64.1403212323220.19717@digraph.polyomino.org.uk |
|---|---|
| State | New |
| Headers | show |
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
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; + } + /* 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; + } /* 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; + 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; + else goto ret; - } - } + } + } ret: return retval; + + ret_huge: + return hhuge * hhuge; + + ret_tiny: + return tiny * tiny; } #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