Fix missing exceptions from exp (bugs 13787, 13922, 14036).

5 files changed, 51 insertions, 52 deletions

diff --git a/sysdeps/i386/fpu/e_expl.S b/sysdeps/i386/fpu/e_expl.S
index a492c29..45c4d07 100644
--- a/sysdeps/i386/fpu/e_expl.S
+++ b/sysdeps/i386/fpu/e_expl.S
@@ -35,6 +35,10 @@ c0:	.byte 0, 0, 0, 0, 0, 0, 0xaa, 0xb8, 0xff, 0x3f
 c1:	.byte 0x20, 0xfa, 0xee, 0xc2, 0x5f, 0x70, 0xa5, 0xec, 0xed, 0x3f
 	.byte 0, 0, 0, 0, 0, 0
 	ASM_SIZE_DIRECTIVE(c1)
+	ASM_TYPE_DIRECTIVE(csat,@object)
+csat:	.byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
+	.byte 0, 0, 0, 0, 0, 0
+	ASM_SIZE_DIRECTIVE(csat)
 
 #ifdef PIC
 # define MO(op) op##@GOTOFF(%ecx)
@@ -53,12 +57,25 @@ ENTRY(__ieee754_expl)
 #ifdef PIC
 	LOAD_PIC_REG (cx)
 #endif
+	movzwl	4+8(%esp), %eax
+	andl	$0x7fff, %eax
+	cmpl	$0x400d, %eax
+	jle	3f
+	/* Overflow, underflow or infinity or NaN as argument.  */
 	fstsw	%ax
 	movb	$0x45, %dh
 	andb	%ah, %dh
 	cmpb	$0x05, %dh
 	je	1f		/* Is +-Inf, jump.    */
-	fldl2e			/* 1  log2(e)         */
+	cmpb	$0x01, %dh
+	je	2f		/* Is +-NaN, jump.    */
+	/* Overflow or underflow; saturate.  */
+	fstp	%st
+	fldt	MO(csat)
+	andb	$2, %ah
+	jz	3f
+	fchs
+3:	fldl2e			/* 1  log2(e)         */
 	fmul	%st(1), %st	/* 1  x log2(e)       */
 	frndint			/* 1  i               */
 	fld	%st(1)		/* 2  x               */
diff --git a/sysdeps/ieee754/dbl-64/w_exp.c b/sysdeps/ieee754/dbl-64/w_exp.c
index aa8ff76..14328a7 100644
--- a/sysdeps/ieee754/dbl-64/w_exp.c
+++ b/sysdeps/ieee754/dbl-64/w_exp.c
@@ -19,27 +19,16 @@
 #include <math.h>
 #include <math_private.h>
 
-static const double
-o_threshold=  7.09782712893383973096e+02,  /* 0x40862E42, 0xFEFA39EF */
-u_threshold= -7.45133219101941108420e+02;  /* 0xc0874910, 0xD52D3051 */
-
-
 /* wrapper exp */
 double
 __exp (double x)
 {
-  if (__builtin_expect (isgreater (x, o_threshold), 0))
-    {
-      if (_LIB_VERSION != _IEEE_)
-	return __kernel_standard (x, x, 6);
-    }
-  else if (__builtin_expect (isless (x, u_threshold), 0))
-    {
-      if (_LIB_VERSION != _IEEE_)
-	return __kernel_standard (x, x, 7);
-    }
+  double z = __ieee754_exp (x);
+  if (__builtin_expect (!__finite (z) || z == 0, 0)
+      && __finite (x) && _LIB_VERSION != _IEEE_)
+    return __kernel_standard (x, x, 6 + !!__signbit (x));
 
-  return __ieee754_exp (x);
+  return z;
 }
 hidden_def (__exp)
 weak_alias (__exp, exp)
diff --git a/sysdeps/ieee754/flt-32/w_expf.c b/sysdeps/ieee754/flt-32/w_expf.c
index bc3b2f6..bfef9e4 100644
--- a/sysdeps/ieee754/flt-32/w_expf.c
+++ b/sysdeps/ieee754/flt-32/w_expf.c
@@ -19,27 +19,16 @@
 #include <math.h>
 #include <math_private.h>
 
-static const float
-o_threshold=  8.8722831726e+01,  /* 0x42b17217 */
-u_threshold= -1.0397208405e+02;  /* 0xc2cff1b5 */
-
-
 /* wrapper expf */
 float
 __expf (float x)
 {
-  if (__builtin_expect (isgreater (x, o_threshold), 0))
-    {
-      if (_LIB_VERSION != _IEEE_)
-	return __kernel_standard_f (x, x, 106);
-    }
-  else if (__builtin_expect (isless (x, u_threshold), 0))
-    {
-      if (_LIB_VERSION != _IEEE_)
-	return __kernel_standard_f (x, x, 107);
-    }
+  float z = __ieee754_expf (x);
+  if (__builtin_expect (!__finitef (z) || z == 0, 0)
+      && __finitef (x) && _LIB_VERSION != _IEEE_)
+    return __kernel_standard_f (x, x, 106 + !!__signbitf (x));
 
-  return __ieee754_expf (x);
+  return z;
 }
 hidden_def (__expf)
 weak_alias (__expf, expf)
diff --git a/sysdeps/ieee754/ldbl-96/w_expl.c b/sysdeps/ieee754/ldbl-96/w_expl.c
index 55c6846..79b10c5 100644
--- a/sysdeps/ieee754/ldbl-96/w_expl.c
+++ b/sysdeps/ieee754/ldbl-96/w_expl.c
@@ -19,29 +19,16 @@
 #include <math.h>
 #include <math_private.h>
 
-static const long double
-o_threshold=  1.135652340629414394949193107797076489134e4,
-  /* 0x400C, 0xB17217F7, 0xD1CF79AC */
-u_threshold= -1.140019167866942050398521670162263001513e4;
-  /* 0x400C, 0xB220C447, 0x69C201E8 */
-
-
 /* wrapper expl */
 long double
 __expl (long double x)
 {
-  if (__builtin_expect (isgreater (x, o_threshold), 0))
-    {
-      if (_LIB_VERSION != _IEEE_)
-	return __kernel_standard_l (x, x, 206);
-    }
-  else if (__builtin_expect (isless (x, u_threshold), 0))
-    {
-      if (_LIB_VERSION != _IEEE_)
-	return __kernel_standard_l (x, x, 207);
-    }
+  long double z = __ieee754_expl (x);
+  if (__builtin_expect (!__finitel (z) || z == 0, 0)
+      && __finitel (x) && _LIB_VERSION != _IEEE_)
+    return __kernel_standard_l (x, x, 206 + !!__signbitl (x));
 
-  return __ieee754_expl (x);
+  return z;
 }
 hidden_def (__expl)
 weak_alias (__expl, expl)
diff --git a/sysdeps/x86_64/fpu/e_expl.S b/sysdeps/x86_64/fpu/e_expl.S
index c2d1d40..d497b28 100644
--- a/sysdeps/x86_64/fpu/e_expl.S
+++ b/sysdeps/x86_64/fpu/e_expl.S
@@ -35,6 +35,10 @@ c0:	.byte 0, 0, 0, 0, 0, 0, 0xaa, 0xb8, 0xff, 0x3f
 c1:	.byte 0x20, 0xfa, 0xee, 0xc2, 0x5f, 0x70, 0xa5, 0xec, 0xed, 0x3f
 	.byte 0, 0, 0, 0, 0, 0
 	ASM_SIZE_DIRECTIVE(c1)
+	ASM_TYPE_DIRECTIVE(csat,@object)
+csat:	.byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
+	.byte 0, 0, 0, 0, 0, 0
+	ASM_SIZE_DIRECTIVE(csat)
 
 #ifdef PIC
 # define MO(op) op##(%rip)
@@ -50,12 +54,25 @@ ENTRY(__ieee754_expl)
    For the i686 the code can be written better.
    -- drepper@cygnus.com.  */
 	fxam			/* Is NaN or +-Inf?  */
+	movzwl	8+8(%rsp), %eax
+	andl	$0x7fff, %eax
+	cmpl	$0x400d, %eax
+	jle	3f
+	/* Overflow, underflow or infinity or NaN as argument.  */
 	fstsw	%ax
 	movb	$0x45, %dh
 	andb	%ah, %dh
 	cmpb	$0x05, %dh
 	je	1f		/* Is +-Inf, jump.    */
-	fldl2e			/* 1  log2(e)         */
+	cmpb	$0x01, %dh
+	je	2f		/* Is +-NaN, jump.    */
+	/* Overflow or underflow; saturate.  */
+	fstp	%st
+	fldt	MO(csat)
+	andb	$2, %ah
+	jz	3f
+	fchs
+3:	fldl2e			/* 1  log2(e)         */
 	fmul	%st(1), %st	/* 1  x log2(e)       */
 	frndint			/* 1  i               */
 	fld	%st(1)		/* 2  x               */