This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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


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

Attachment: pgpeO6e2NuPdH.pgp
Description: PGP signature


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]