Similar to various other bugs in this area, some asin 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, powerpc and mips64. Committed.
(auto-libm-test-out diffs omitted below.)
2015-02-26 Joseph Myers <joseph@codesourcery.com>
[BZ #16351]
* sysdeps/i386/fpu/e_asin.S (dbl_min): New object.
(MO): New macro.
(__ieee754_asin): Force underflow exception for results with small
absolute value.
* sysdeps/i386/fpu/e_asinf.S (flt_min): New object.
(MO): New macro.
(__ieee754_asinf): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/dbl-64/e_asin.c: Include <float.h> and <math.h>.
(__ieee754_asin): Force underflow exception for results with small
absolute value.
* sysdeps/ieee754/flt-32/e_asinf.c: Include <float.h>.
(__ieee754_asinf): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/ldbl-128/e_asinl.c: Include <float.h>.
(__ieee754_asinl): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/ldbl-128ibm/e_asinl.c: Include <float.h>.
(__ieee754_asinl): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/ldbl-96/e_asinl.c: Include <float.h>.
(__ieee754_asinl): Force underflow exception for results with
small absolute value.
* sysdeps/x86_64/fpu/multiarch/e_asin.c [HAVE_FMA4_SUPPORT]:
Include <math.h>.
* math/auto-libm-test-in: Do not mark underflow exceptions as
possibly missing for bug 16351.
* math/auto-libm-test-out: Regenerated.
@@ -91,11 +91,10 @@ asin 0x0.ffffffffffffp0
asin -0x0.ffffffffffffp0
asin 0x0.ffffffffffffffffp0
asin -0x0.ffffffffffffffffp0
-# Bug 16351: underflow exception may be missing.
-asin min missing-underflow
-asin -min missing-underflow
-asin min_subnorm missing-underflow
-asin -min_subnorm missing-underflow
+asin min
+asin -min
+asin min_subnorm
+asin -min_subnorm
asinh 0
asinh -0
@@ -7,8 +7,26 @@
RCSID("$NetBSD: e_asin.S,v 1.4 1995/05/08 23:45:40 jtc Exp $")
+ .section .rodata.cst8,"aM",@progbits,8
+
+ .p2align 3
+ .type dbl_min,@object
+dbl_min: .byte 0, 0, 0, 0, 0, 0, 0x10, 0
+ ASM_SIZE_DIRECTIVE(dbl_min)
+
+#ifdef PIC
+# define MO(op) op##@GOTOFF(%ecx)
+#else
+# define MO(op) op
+#endif
+
+ .text
+
/* asin = atan (x / sqrt((1-x) (1+x))) */
ENTRY(__ieee754_asin)
+#ifdef PIC
+ LOAD_PIC_REG (cx)
+#endif
fldl 4(%esp) /* x */
fld %st
fld1 /* 1 : x : x */
@@ -18,6 +36,20 @@ ENTRY(__ieee754_asin)
fmulp /* 1 - x^2 */
fsqrt /* sqrt (1 - x^2) */
fpatan
- ret
+ fldl MO(dbl_min)
+ fld %st(1)
+ fabs
+ fucompp
+ fnstsw
+ sahf
+ jnc 1f
+ subl $8, %esp
+ cfi_adjust_cfa_offset (8)
+ fld %st(0)
+ fmul %st(0)
+ fstpl (%esp)
+ addl $8, %esp
+ cfi_adjust_cfa_offset (-8)
+1: ret
END (__ieee754_asin)
strong_alias (__ieee754_asin, __asin_finite)
@@ -8,8 +8,26 @@
RCSID("$NetBSD: $")
+ .section .rodata.cst4,"aM",@progbits,4
+
+ .p2align 2
+ .type flt_min,@object
+flt_min: .byte 0, 0, 0x80, 0
+ ASM_SIZE_DIRECTIVE(flt_min)
+
+#ifdef PIC
+# define MO(op) op##@GOTOFF(%ecx)
+#else
+# define MO(op) op
+#endif
+
+ .text
+
/* asin = atan (x / sqrt(1 - x^2)) */
ENTRY(__ieee754_asinf)
+#ifdef PIC
+ LOAD_PIC_REG (cx)
+#endif
flds 4(%esp) /* x */
fld %st
fmul %st(0) /* x^2 */
@@ -17,6 +35,20 @@ ENTRY(__ieee754_asinf)
fsubp /* 1 - x^2 */
fsqrt /* sqrt (1 - x^2) */
fpatan
- ret
+ flds MO(flt_min)
+ fld %st(1)
+ fabs
+ fucompp
+ fnstsw
+ sahf
+ jnc 1f
+ subl $4, %esp
+ cfi_adjust_cfa_offset (4)
+ fld %st(0)
+ fmul %st(0)
+ fstps (%esp)
+ addl $4, %esp
+ cfi_adjust_cfa_offset (-4)
+1: ret
END (__ieee754_asinf)
strong_alias (__ieee754_asinf, __asinf_finite)
@@ -39,6 +39,8 @@
#include "powtwo.tbl"
#include "MathLib.h"
#include "uasncs.h"
+#include <float.h>
+#include <math.h>
#include <math_private.h>
#ifndef SECTION
@@ -67,7 +69,15 @@ __ieee754_asin(double x){
m = u.i[HIGH_HALF];
k = 0x7fffffff&m; /* no sign */
- if (k < 0x3e500000) return x; /* for x->0 => sin(x)=x */
+ if (k < 0x3e500000)
+ {
+ if (fabs (x) < DBL_MIN)
+ {
+ double force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
+ return x; /* for x->0 => sin(x)=x */
+ }
/*----------------------2^-26 <= |x| < 2^ -3 -----------------*/
else
if (k < 0x3fc00000) {
@@ -39,6 +39,7 @@
static char rcsid[] = "$NetBSD: e_asinf.c,v 1.5 1995/05/12 04:57:25 jtc Exp $";
#endif
+#include <float.h>
#include <math.h>
#include <math_private.h>
@@ -72,6 +73,11 @@ float __ieee754_asinf(float x)
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix<0x3f000000) { /* |x|<0.5 */
if(ix<0x32000000) { /* if |x| < 2**-27 */
+ if (fabsf (x) < FLT_MIN)
+ {
+ float force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
if(huge+x>one) return x;/* return x with inexact if x!=0*/
} else {
t = x*x;
@@ -59,6 +59,7 @@
*/
+#include <float.h>
#include <math.h>
#include <math_private.h>
long double sqrtl (long double);
@@ -152,6 +153,11 @@ __ieee754_asinl (long double x)
{
if (ix < 0x3fc60000) /* |x| < 2**-57 */
{
+ if (fabsl (x) < LDBL_MIN)
+ {
+ long double force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
if (huge + x > one)
return x; /* return x with inexact if x!=0 */
}
@@ -59,6 +59,7 @@
*/
+#include <float.h>
#include <math.h>
#include <math_private.h>
long double sqrtl (long double);
@@ -146,6 +147,11 @@ __ieee754_asinl (long double x)
{
if (a < 6.938893903907228e-18L) /* |x| < 2**-57 */
{
+ if (fabsl (x) < LDBL_MIN)
+ {
+ long double force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
if (huge + x > one)
return x; /* return x with inexact if x!=0 */
}
@@ -58,6 +58,7 @@
*/
+#include <float.h>
#include <math.h>
#include <math_private.h>
@@ -111,6 +112,11 @@ __ieee754_asinl (long double x)
{ /* |x|<0.5 */
if (ix < 0x3fde8000)
{ /* if |x| < 2**-33 */
+ if (fabsl (x) < LDBL_MIN)
+ {
+ long double force_underflow = x * x;
+ math_force_eval (force_underflow);
+ }
if (huge + x > one)
return x; /* return x with inexact if x!=0 */
}
@@ -1,5 +1,6 @@
#ifdef HAVE_FMA4_SUPPORT
# include <init-arch.h>
+# include <math.h>
# include <math_private.h>
extern double __ieee754_acos_sse2 (double);