Fix dbl-64 exp overflow/underflow in non-default rounding modes (bug 16284)
diff mbox

Message ID Pine.LNX.4.64.1403212323220.19717@digraph.polyomino.org.uk
State New
Headers show

Commit Message

Joseph Myers March 21, 2014, 11:24 p.m. UTC
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.

Comments

Siddhesh Poyarekar March 24, 2014, 9:16 a.m. UTC | #1
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

Patch
diff mbox

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