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

8 files changed, 87 insertions, 64 deletions

diff --git a/ChangeLog b/ChangeLog
index 5ae8e58404..06fe04f328 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,28 @@
 2012-05-05  Joseph Myers  <joseph@codesourcery.com>
 
+	[BZ #13787]
+	[BZ #13922]
+	[BZ #14036]
+	* sysdeps/i386/fpu/e_expl.S (csat): New constant.
+	(__ieee754_expl): Allow for and saturate large arguments.
+	* sysdeps/ieee754/dbl-64/w_exp.c (o_threshold): Remove variable.
+	(u_threshold): Likewise.
+	(__exp): Call __ieee754_exp before checking for overflow and
+	underflow.
+	* sysdeps/ieee754/flt-32/w_expf.c (o_threshold): Remove variable.
+	(u_threshold): Likewise.
+	(__expf): Call __ieee754_expf before checking for overflow and
+	underflow.
+	* sysdeps/ieee754/ldbl-96/w_expl.c (o_threshold): Remove variable.
+	(u_threshold): Likewise.
+	(__expl): Call __ieee754_expl before checking for overflow and
+	underflow.
+	* sysdeps/x86_64/fpu/e_expl.S (csat): New constant.
+	(__ieee754_expl): Allow for and saturate large arguments.
+	* math/libm-test.inc (exp_test): Add another test.  Do not allow
+	missing overflow exception on overflow.
+	(expm1_test): Do not allow missing overflow exception on overflow.
+
 	* sysdeps/i386/fpu/e_expl.c: Move to ...
 	* sysdeps/i386/fpu/e_expl.S: ... here.  Write directly in assembly
 	rather than using inline asm.
diff --git a/NEWS b/NEWS
index 0873b2f03d..f75fed6a44 100644
--- a/NEWS
+++ b/NEWS
@@ -18,12 +18,12 @@ Version 2.16
   12047, 12340, 13058, 13525, 13526, 13527, 13528, 13529, 13530, 13531,
   13532, 13533, 13547, 13551, 13552, 13553, 13555, 13559, 13566, 13583,
   13592, 13618, 13637, 13656, 13658, 13673, 13691, 13695, 13704, 13705,
-  13706, 13726, 13738, 13739, 13758, 13760, 13761, 13775, 13786, 13792,
-  13806, 13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871,
-  13872, 13873, 13879, 13883, 13886, 13892, 13895, 13908, 13910, 13911,
-  13912, 13913, 13915, 13916, 13917, 13918, 13919, 13920, 13921, 13924,
-  13926, 13927, 13928, 13938, 13941, 13942, 13963, 13967, 13970, 13973,
-  14027, 14033, 14034, 14040, 14049, 14055
+  13706, 13726, 13738, 13739, 13758, 13760, 13761, 13775, 13786, 13787,
+  13792, 13806, 13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854,
+  13871, 13872, 13873, 13879, 13883, 13886, 13892, 13895, 13908, 13910,
+  13911, 13912, 13913, 13915, 13916, 13917, 13918, 13919, 13920, 13921,
+  13922, 13924, 13926, 13927, 13928, 13938, 13941, 13942, 13963, 13967,
+  13970, 13973, 14027, 14033, 14034, 14040, 14049, 14055
 
 * ISO C11 support:
 
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 0875e2c467..fce27f31d2 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -3325,8 +3325,11 @@ exp_test (void)
   TEST_f_f (exp, 1000.0L, 0.197007111401704699388887935224332313e435L);
 #endif
 
-  /* Bug 13922: OVERFLOW exception may be missing.  */
-  TEST_f_f (exp, max_value, plus_infty, OVERFLOW_EXCEPTION_OK);
+#if !(defined TEST_LDOUBLE && LDBL_MAX_EXP > 1024)
+  TEST_f_f (exp, 710, plus_infty, OVERFLOW_EXCEPTION);
+#endif
+  TEST_f_f (exp, 1e5, plus_infty, OVERFLOW_EXCEPTION);
+  TEST_f_f (exp, max_value, plus_infty, OVERFLOW_EXCEPTION);
   TEST_f_f (exp, -max_value, 0);
 
   END (exp);
@@ -3547,11 +3550,9 @@ expm1_test (void)
 #endif
 
   errno = 0;
-  /* Bug 13787: OVERFLOW exception may be missing.  */
-  TEST_f_f (expm1, 100000.0, plus_infty, OVERFLOW_EXCEPTION_OK);
+  TEST_f_f (expm1, 100000.0, plus_infty, OVERFLOW_EXCEPTION);
   check_int ("errno for expm1(large) == ERANGE", errno, ERANGE, 0, 0, 0);
-  /* Bug 13787: OVERFLOW exception may be missing.  */
-  TEST_f_f (expm1, max_value, plus_infty, OVERFLOW_EXCEPTION_OK);
+  TEST_f_f (expm1, max_value, plus_infty, OVERFLOW_EXCEPTION);
 #ifndef TEST_LDOUBLE /* Bug 13923.  */
   TEST_f_f (expm1, -max_value, -1);
 #endif
diff --git a/sysdeps/i386/fpu/e_expl.S b/sysdeps/i386/fpu/e_expl.S
index a492c29a75..45c4d07539 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 aa8ff7689f..14328a7b4b 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 bc3b2f6790..bfef9e4d24 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 55c68462bd..79b10c5756 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 c2d1d40ba5..d497b28973 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               */