Fix sinh missing underflows (bug 16519).

Similar to various other bugs in this area, some sinh implementations do not raise the underflow exception for subnormal arguments, when the result is tiny and inexact. This patch forces the exception in a similar way to previous fixes. Tested for x86_64, x86, mips64 and powerpc. [BZ #16519] * sysdeps/ieee754/dbl-64/e_sinh.c: Include <float.h>. (__ieee754_sinh): Force underflow exception for arguments with small absolute value. * sysdeps/ieee754/flt-32/e_sinhf.c: Include <float.h>. (__ieee754_sinhf): Force underflow exception for arguments with small absolute value. * sysdeps/ieee754/ldbl-128/e_sinhl.c: Include <float.h>. (__ieee754_sinhl): Force underflow exception for arguments with small absolute value. * sysdeps/ieee754/ldbl-128ibm/e_sinhl.c: Include <float.h>. (__ieee754_sinhl): Force underflow exception for arguments with small absolute value. * sysdeps/ieee754/ldbl-96/e_sinhl.c: Include <float.h>. (__ieee754_sinhl): Force underflow exception for arguments with small absolute value. * math/auto-libm-test-in: Add more tests of sinh. * math/auto-libm-test-out: Regenerated. * sysdeps/i386/fpu/libm-test-ulps: Update.
git-mirror · Aug 6, 2015 · 5e29dd5 · 5e29dd5
1 parent 2ba30a1
commit 5e29dd5
Show file tree

Hide file tree

Showing 10 changed files with 392 additions and 8 deletions.
diff --git a/ChangeLog b/ChangeLog
@@ -1,3 +1,25 @@
+2015-08-06  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #16519]
+	* sysdeps/ieee754/dbl-64/e_sinh.c: Include <float.h>.
+	(__ieee754_sinh): Force underflow exception for arguments with
+	small absolute value.
+	* sysdeps/ieee754/flt-32/e_sinhf.c: Include <float.h>.
+	(__ieee754_sinhf): Force underflow exception for arguments with
+	small absolute value.
+	* sysdeps/ieee754/ldbl-128/e_sinhl.c: Include <float.h>.
+	(__ieee754_sinhl): Force underflow exception for arguments with
+	small absolute value.
+	* sysdeps/ieee754/ldbl-128ibm/e_sinhl.c: Include <float.h>.
+	(__ieee754_sinhl): Force underflow exception for arguments with
+	small absolute value.
+	* sysdeps/ieee754/ldbl-96/e_sinhl.c: Include <float.h>.
+	(__ieee754_sinhl): Force underflow exception for arguments with
+	small absolute value.
+	* math/auto-libm-test-in: Add more tests of sinh.
+	* math/auto-libm-test-out: Regenerated.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+
 2015-08-06  Paul Pluzhnikov  <ppluzhnikov@google.com>
 
 	* libio/oldfileops.c: Include unistd.h.

diff --git a/NEWS b/NEWS
@@ -8,7 +8,7 @@ using `glibc' in the "product" field.
 Version 2.23
 
 * The following bugs are resolved with this release:
-  18265, 18525, 18647, 18661.
+  16519, 18265, 18525, 18647, 18661.
 
 Version 2.22
 

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
@@ -2528,6 +2528,10 @@ sinh 0x7.6e259d2436fc4p-4
 sinh 0x3.d6e088p-4
 sinh -0x7.688eap-4
 sinh -0xd.dce79p-4
+sinh min
+sinh -min
+sinh min_subnorm
+sinh -min_subnorm
 
 sqrt 0
 sqrt -0

diff --git a/math/auto-libm-test-out b/math/auto-libm-test-out
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
@@ -1864,6 +1864,10 @@ ildouble: 2
 ldouble: 2
 
 Function: "sinh_downward":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
 ildouble: 4
 ldouble: 5
 
@@ -1874,6 +1878,8 @@ ldouble: 4
 Function: "sinh_upward":
 double: 1
 float: 1
+idouble: 1
+ifloat: 1
 ildouble: 4
 ldouble: 5
 

diff --git a/sysdeps/ieee754/dbl-64/e_sinh.c b/sysdeps/ieee754/dbl-64/e_sinh.c
@@ -32,6 +32,7 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $";
  *	only sinh(0)=0 is exact for finite x.
  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -58,10 +59,16 @@ __ieee754_sinh (double x)
   /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
   if (ix < 0x40360000)                  /* |x|<22 */
     {
-      if (__glibc_unlikely (ix < 0x3e300000))              /* |x|<2**-28 */
+      if (__glibc_unlikely (ix < 0x3e300000)) {            /* |x|<2**-28 */
+	if (fabs (x) < DBL_MIN)
+	  {
+	    double force_underflow = x * x;
+	    math_force_eval (force_underflow);
+	  }
 	if (shuge + x > one)
 	  return x;
-      /* sinh(tiny) = tiny with inexact */
+	/* sinh(tiny) = tiny with inexact */
+      }
       t = __expm1 (fabs (x));
       if (ix < 0x3ff00000)
 	return h * (2.0 * t - t * t / (t + one));

diff --git a/sysdeps/ieee754/flt-32/e_sinhf.c b/sysdeps/ieee754/flt-32/e_sinhf.c
@@ -13,6 +13,7 @@
  * ====================================================
  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -34,8 +35,14 @@ __ieee754_sinhf(float x)
 	if (jx<0) h = -h;
     /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
 	if (ix < 0x41b00000) {		/* |x|<22 */
-	    if (__builtin_expect(ix<0x31800000, 0))	/* |x|<2**-28 */
+	    if (__builtin_expect(ix<0x31800000, 0)) {	/* |x|<2**-28 */
+		if (fabsf (x) < FLT_MIN)
+		  {
+		    float force_underflow = x * x;
+		    math_force_eval (force_underflow);
+		  }
 		if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
+	    }
 	    t = __expm1f(fabsf(x));
 	    if(ix<0x3f800000) return h*((float)2.0*t-t*t/(t+one));
 	    return h*(t+t/(t+one));

diff --git a/sysdeps/ieee754/ldbl-128/e_sinhl.c b/sysdeps/ieee754/ldbl-128/e_sinhl.c
@@ -53,6 +53,7 @@
  *      only sinhl(0)=0 is exact for finite x.
  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -86,8 +87,15 @@ __ieee754_sinhl (long double x)
   if (ix <= 0x40044000)
     {
       if (ix < 0x3fc60000) /* |x| < 2^-57 */
-	if (shuge + x > one)
-	  return x;		/* sinh(tiny) = tiny with inexact */
+	{
+	  if (fabsl (x) < LDBL_MIN)
+	    {
+	      long double force_underflow = x * x;
+	      math_force_eval (force_underflow);
+	    }
+	  if (shuge + x > one)
+	    return x;		/* sinh(tiny) = tiny with inexact */
+	}
       t = __expm1l (u.value);
       if (ix < 0x3fff0000)
 	return h * (2.0 * t - t * t / (t + one));

diff --git a/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c b/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
@@ -28,6 +28,7 @@
  *	only sinh(0)=0 is exact for finite x.
  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -52,8 +53,14 @@ __ieee754_sinhl(long double x)
 	if (jx<0) h = -h;
     /* |x| in [0,40], return sign(x)*0.5*(E+E/(E+1))) */
 	if (ix < 0x4044000000000000LL) {	/* |x|<40 */
-	    if (ix<0x3e20000000000000LL)	/* |x|<2**-29 */
+	    if (ix<0x3e20000000000000LL) {	/* |x|<2**-29 */
+		if (fabsl (x) < LDBL_MIN)
+		  {
+		    long double force_underflow = x * x;
+		    math_force_eval (force_underflow);
+		  }
 		if(shuge+x>one) return x;/* sinhl(tiny) = tiny with inexact */
+	    }
 	    t = __expm1l(fabsl(x));
 	    if(ix<0x3ff0000000000000LL) return h*(2.0*t-t*t/(t+one));
 	    w = t/(t+one);

diff --git a/sysdeps/ieee754/ldbl-96/e_sinhl.c b/sysdeps/ieee754/ldbl-96/e_sinhl.c
@@ -36,6 +36,7 @@ static char rcsid[] = "$NetBSD: $";
  *	only sinhl(0)=0 is exact for finite x.
  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -58,8 +59,14 @@ __ieee754_sinhl(long double x)
 	if (jx & 0x8000) h = -h;
     /* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */
 	if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x|<25 */
-	    if (ix<0x3fdf)		/* |x|<2**-32 */
+	    if (ix<0x3fdf) {		/* |x|<2**-32 */
+		if (fabsl (x) < LDBL_MIN)
+		  {
+		    long double force_underflow = x * x;
+		    math_force_eval (force_underflow);
+		  }
 		if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
+	    }
 	    t = __expm1l(fabsl(x));
 	    if(ix<0x3fff) return h*(2.0*t-t*t/(t+one));
 	    return h*(t+t/(t+one));