This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix sin, sincos missing underflows (bug 16526, bug 16538) [committed]
- From: Joseph Myers <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Tue, 23 Jun 2015 22:25:41 +0000
- Subject: Fix sin, sincos missing underflows (bug 16526, bug 16538) [committed]
- Authentication-results: sourceware.org; auth=none
Similar to various other bugs in this area, some sin and sincos
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. Committed.
(auto-libm-test-out diffs omitted below.)
2015-06-23 Joseph Myers <joseph@codesourcery.com>
[BZ #16526]
[BZ #16538]
* sysdeps/ieee754/dbl-64/s_sin.c: Include <float.h>.
(__sin): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/flt-32/k_sinf.c: Include <float.h>.
(__kernel_sinf): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128/k_sincosl.c: Include <float.h>.
(__kernel_sincosl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128/k_sinl.c: Include <float.h>.
(__kernel_sinl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128ibm/k_sincosl.c: Include <float.h>.
(__kernel_sincosl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128ibm/k_sinl.c: Include <float.h>.
(__kernel_sinl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-96/k_sinl.c: Include <float.h>.
(__kernel_sinl): Force underflow exception for arguments with
small absolute value.
* sysdeps/powerpc/fpu/k_sinf.c: Include <float.h>.
(__kernel_sinf): Force underflow exception for arguments with
small absolute value.
* math/auto-libm-test-in: Add more tests of sin and sincos.
* math/auto-libm-test-out: Regenerated.
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index d1502ad..6545583 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -2317,6 +2317,10 @@ sin 0x1.2001469775ce6p32
sin -0x3.3de320f6be87ep+1020
sin 0xe.9f1e5bc3bb88p+112
sin 0x4.7857dp+68
+sin min
+sin -min
+sin min_subnorm
+sin -min_subnorm
sincos 0
sincos -0
@@ -2339,6 +2343,10 @@ sincos 0x1p+28
sincos -0x3.3de320f6be87ep+1020
sincos 0xe.9f1e5bc3bb88p+112
sincos 0x4.7857dp+68
+sincos min
+sincos -min
+sincos min_subnorm
+sincos -min_subnorm
sinh 0
sinh -0
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index ea89ad5..eff120e 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -48,6 +48,7 @@
#include <errno.h>
+#include <float.h>
#include "endian.h"
#include "mydefs.h"
#include "usncs.h"
@@ -295,7 +296,14 @@ __sin (double x)
m = u.i[HIGH_HALF];
k = 0x7fffffff & m; /* no sign */
if (k < 0x3e500000) /* if x->0 =>sin(x)=x */
- retval = x;
+ {
+ if (fabs (x) < DBL_MIN)
+ {
+ double force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
+ retval = x;
+ }
/*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
else if (k < 0x3fd00000)
{
diff --git a/sysdeps/ieee754/flt-32/k_sinf.c b/sysdeps/ieee754/flt-32/k_sinf.c
index 0bafd83..0c98a2a 100644
--- a/sysdeps/ieee754/flt-32/k_sinf.c
+++ b/sysdeps/ieee754/flt-32/k_sinf.c
@@ -17,6 +17,7 @@
static char rcsid[] = "$NetBSD: k_sinf.c,v 1.4 1995/05/10 20:46:33 jtc Exp $";
#endif
+#include <float.h>
#include <math.h>
#include <math_private.h>
@@ -36,7 +37,15 @@ float __kernel_sinf(float x, float y, int iy)
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff; /* high word of x */
if(ix<0x32000000) /* |x| < 2**-27 */
- {if((int)x==0) return x;} /* generate inexact */
+ {
+ if (fabsf (x) < FLT_MIN)
+ {
+ float force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
+ if ((int) x == 0)
+ return x; /* generate inexact */
+ }
z = x*x;
v = z*x;
r = S2+z*(S3+z*(S4+z*(S5+z*S6)));
diff --git a/sysdeps/ieee754/ldbl-128/k_sincosl.c b/sysdeps/ieee754/ldbl-128/k_sincosl.c
index 2f66dee..7b5c4b0 100644
--- a/sysdeps/ieee754/ldbl-128/k_sincosl.c
+++ b/sysdeps/ieee754/ldbl-128/k_sincosl.c
@@ -17,6 +17,7 @@
License along with the GNU C Library; if not, see
<http://www.gnu.org/licenses/>. */
+#include <float.h>
#include <math.h>
#include <math_private.h>
@@ -109,12 +110,19 @@ __kernel_sincosl(long double x, long double y, long double *sinx, long double *c
/* Argument is small enough to approximate it by a Chebyshev
polynomial of degree 16(17). */
if (tix < 0x3fc60000) /* |x| < 2^-57 */
- if (!((int)x)) /* generate inexact */
- {
- *sinx = x;
- *cosx = ONE;
- return;
- }
+ {
+ if (fabsl (x) < LDBL_MIN)
+ {
+ long double force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
+ if (!((int)x)) /* generate inexact */
+ {
+ *sinx = x;
+ *cosx = ONE;
+ return;
+ }
+ }
z = x * x;
*sinx = x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+
z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8)))))))));
diff --git a/sysdeps/ieee754/ldbl-128/k_sinl.c b/sysdeps/ieee754/ldbl-128/k_sinl.c
index b15521b..04d539f 100644
--- a/sysdeps/ieee754/ldbl-128/k_sinl.c
+++ b/sysdeps/ieee754/ldbl-128/k_sinl.c
@@ -17,6 +17,7 @@
License along with the GNU C Library; if not, see
<http://www.gnu.org/licenses/>. */
+#include <float.h>
#include <math.h>
#include <math_private.h>
@@ -90,7 +91,14 @@ __kernel_sinl(long double x, long double y, int iy)
/* Argument is small enough to approximate it by a Chebyshev
polynomial of degree 17. */
if (tix < 0x3fc60000) /* |x| < 2^-57 */
- if (!((int)x)) return x; /* generate inexact */
+ {
+ if (fabsl (x) < LDBL_MIN)
+ {
+ long double force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
+ if (!((int)x)) return x; /* generate inexact */
+ }
z = x * x;
return x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+
z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8)))))))));
diff --git a/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c b/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c
index f831e3f..0a76e1c 100644
--- a/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c
@@ -17,6 +17,7 @@
License along with the GNU C Library; if not, see
<http://www.gnu.org/licenses/>. */
+#include <float.h>
#include <math.h>
#include <math_private.h>
@@ -112,12 +113,19 @@ __kernel_sincosl(long double x, long double y, long double *sinx, long double *c
/* Argument is small enough to approximate it by a Chebyshev
polynomial of degree 16(17). */
if (tix < 0x3c600000) /* |x| < 2^-57 */
- if (!((int)x)) /* generate inexact */
- {
- *sinx = x;
- *cosx = ONE;
- return;
- }
+ {
+ if (fabsl (x) < LDBL_MIN)
+ {
+ long double force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
+ if (!((int)x)) /* generate inexact */
+ {
+ *sinx = x;
+ *cosx = ONE;
+ return;
+ }
+ }
z = x * x;
*sinx = x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+
z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8)))))))));
diff --git a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
index d6602fe..2050cd2 100644
--- a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
@@ -17,6 +17,7 @@
License along with the GNU C Library; if not, see
<http://www.gnu.org/licenses/>. */
+#include <float.h>
#include <math.h>
#include <math_private.h>
@@ -93,7 +94,14 @@ __kernel_sinl(long double x, long double y, int iy)
/* Argument is small enough to approximate it by a Chebyshev
polynomial of degree 17. */
if (tix < 0x3c600000) /* |x| < 2^-57 */
- if (!((int)x)) return x; /* generate inexact */
+ {
+ if (fabsl (x) < LDBL_MIN)
+ {
+ long double force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
+ if (!((int)x)) return x; /* generate inexact */
+ }
z = x * x;
return x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+
z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8)))))))));
diff --git a/sysdeps/ieee754/ldbl-96/k_sinl.c b/sysdeps/ieee754/ldbl-96/k_sinl.c
index 179262a..b7b5ae3 100644
--- a/sysdeps/ieee754/ldbl-96/k_sinl.c
+++ b/sysdeps/ieee754/ldbl-96/k_sinl.c
@@ -20,6 +20,7 @@
/* The polynomials have not been optimized for extended-precision and
may contain more terms than needed. */
+#include <float.h>
#include <math.h>
#include <math_private.h>
@@ -94,7 +95,14 @@ __kernel_sinl(long double x, long double y, int iy)
/* Argument is small enough to approximate it by a Chebyshev
polynomial of degree 17. */
if (absx < 0x1p-33L)
- if (!((int)x)) return x; /* generate inexact */
+ {
+ if (fabsl (x) < LDBL_MIN)
+ {
+ long double force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
+ if (!((int)x)) return x; /* generate inexact */
+ }
z = x * x;
return x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+
z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8)))))))));
diff --git a/sysdeps/powerpc/fpu/k_sinf.c b/sysdeps/powerpc/fpu/k_sinf.c
index 64df0c9..e2850df 100644
--- a/sysdeps/powerpc/fpu/k_sinf.c
+++ b/sysdeps/powerpc/fpu/k_sinf.c
@@ -17,6 +17,7 @@
License along with the GNU C Library; see the file COPYING.LIB. If
not, see <http://www.gnu.org/licenses/>. */
+#include <float.h>
#include <math.h>
#include <fenv.h>
#include <math_private.h>
@@ -40,7 +41,10 @@ __kernel_sinf (float x, float y, int iy)
ix = __builtin_fabsf (x);
if (ix < twom27)
{ /* |x| < 2**-27 */
- __feraiseexcept (FE_INEXACT);
+ if (ix < FLT_MIN && ix != 0.0f)
+ __feraiseexcept (FE_UNDERFLOW|FE_INEXACT);
+ else
+ __feraiseexcept (FE_INEXACT);
return x;
}
z = x * x;
--
Joseph S. Myers
joseph@codesourcery.com