Fix asin missing underflows (bug 16351).

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.

	[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.

8 files changed, 102 insertions, 3 deletions

diff --git a/sysdeps/i386/fpu/e_asin.S b/sysdeps/i386/fpu/e_asin.S
index 9a44cb6..ea8cb6f 100644
--- a/sysdeps/i386/fpu/e_asin.S
+++ b/sysdeps/i386/fpu/e_asin.S
@@ -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)
diff --git a/sysdeps/i386/fpu/e_asinf.S b/sysdeps/i386/fpu/e_asinf.S
index 5c1065d..f7bda77 100644
--- a/sysdeps/i386/fpu/e_asinf.S
+++ b/sysdeps/i386/fpu/e_asinf.S
@@ -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)
diff --git a/sysdeps/ieee754/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c
index e90f47c..b87e557 100644
--- a/sysdeps/ieee754/dbl-64/e_asin.c
+++ b/sysdeps/ieee754/dbl-64/e_asin.c
@@ -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) {
diff --git a/sysdeps/ieee754/flt-32/e_asinf.c b/sysdeps/ieee754/flt-32/e_asinf.c
index 91ca9af..00bad42 100644
--- a/sysdeps/ieee754/flt-32/e_asinf.c
+++ b/sysdeps/ieee754/flt-32/e_asinf.c
@@ -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;
diff --git a/sysdeps/ieee754/ldbl-128/e_asinl.c b/sysdeps/ieee754/ldbl-128/e_asinl.c
index cd1a8da..353603d 100644
--- a/sysdeps/ieee754/ldbl-128/e_asinl.c
+++ b/sysdeps/ieee754/ldbl-128/e_asinl.c
@@ -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 */
 	}
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_asinl.c b/sysdeps/ieee754/ldbl-128ibm/e_asinl.c
index dece118..00386db 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_asinl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_asinl.c
@@ -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 */
 	}
diff --git a/sysdeps/ieee754/ldbl-96/e_asinl.c b/sysdeps/ieee754/ldbl-96/e_asinl.c
index c1ffa3e..2973bf0 100644
--- a/sysdeps/ieee754/ldbl-96/e_asinl.c
+++ b/sysdeps/ieee754/ldbl-96/e_asinl.c
@@ -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 */
 	}
diff --git a/sysdeps/x86_64/fpu/multiarch/e_asin.c b/sysdeps/x86_64/fpu/multiarch/e_asin.c
index e742a9c..55865c0 100644
--- a/sysdeps/x86_64/fpu/multiarch/e_asin.c
+++ b/sysdeps/x86_64/fpu/multiarch/e_asin.c
@@ -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);