Fix cexp overflow (bug 13892).

4 files changed, 98 insertions, 24 deletions

diff --git a/math/libm-test.inc b/math/libm-test.inc
index fad767d..3851855 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -1909,6 +1909,35 @@ cexp_test (void)
   TEST_c_c (cexp, -10000, 0x1p16383L, 1.045876464564882298442774542991176546722e-4343L, 4.421154026488516836023811173959413420548e-4344L);
 #endif
 
+  TEST_c_c (cexp, 88.75, 0.75, 2.558360358486542817001900410314204322891e38L, 2.383359453227311447654736314679677655100e38L);
+  TEST_c_c (cexp, -95, 0.75, 4.039714446238306526889476684000081624047e-42L, 3.763383677300535390271646960780570275931e-42L);
+
+#ifndef TEST_FLOAT
+  TEST_c_c (cexp, 709.8125, 0.75, 1.355121963080879535248452862759108365762e308L, 1.262426823598609432507811340856186873507e308L);
+  TEST_c_c (cexp, -720, 0.75, 1.486960657116368433685753325516638551722e-313L, 1.385247284245720590980701226843815229385e-313L);
+#endif
+
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (cexp, 11356.5625, 0.75, 9.052188470850960144814815984311663764287e4931L, 8.432986734191301036267148978260970230200e4931L);
+  TEST_c_c (cexp, -11370, 0.75, 8.631121063182211587489310508568170739592e-4939L, 8.040721827809267291427062346918413482824e-4939L);
+#endif
+
+#ifdef TEST_FLOAT
+  TEST_c_c (cexp, 180, 0x1p-149, plus_infty, 2.087071793345235105931967606907855310664e33L, OVERFLOW_EXCEPTION);
+#endif
+
+#if defined TEST_DOUBLE || (defined TEST_LDOUBLE && LDBL_MAX_EXP == 1024)
+  TEST_c_c (cexp, 1440, 0x1p-1074, plus_infty, 1.196295853897226111293303155636183216483e302L, OVERFLOW_EXCEPTION);
+#endif
+
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (cexp, 22730, 0x1p-16434L, plus_infty, 2.435706297811211974162115164702304105374e4924L, OVERFLOW_EXCEPTION);
+#endif
+
+  TEST_c_c (cexp, 1e6, 0, plus_infty, 0, OVERFLOW_EXCEPTION);
+  TEST_c_c (cexp, 1e6, min_value, plus_infty, plus_infty, OVERFLOW_EXCEPTION);
+  TEST_c_c (cexp, 1e6, -min_value, plus_infty, minus_infty, OVERFLOW_EXCEPTION);
+
   END (cexp, complex);
 }
 
diff --git a/math/s_cexp.c b/math/s_cexp.c
index 82fe814..1d7a5a2 100644
--- a/math/s_cexp.c
+++ b/math/s_cexp.c
@@ -1,5 +1,5 @@
 /* Return value of complex exponential function for double complex value.
-   Copyright (C) 1997, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -21,7 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ double
 __cexp (__complex__ double x)
@@ -36,20 +36,35 @@ __cexp (__complex__ double x)
       if (__builtin_expect (icls >= FP_ZERO, 1))
 	{
 	  /* Imaginary part is finite.  */
-	  double exp_val = __ieee754_exp (__real__ x);
+	  const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2);
 	  double sinix, cosix;
 
 	  __sincos (__imag__ x, &sinix, &cosix);
 
-	  if (isfinite (exp_val))
+	  if (__real__ x > t)
 	    {
-	      __real__ retval = exp_val * cosix;
-	      __imag__ retval = exp_val * sinix;
+	      double exp_t = __ieee754_exp (t);
+	      __real__ x -= t;
+	      sinix *= exp_t;
+	      cosix *= exp_t;
+	      if (__real__ x > t)
+		{
+		  __real__ x -= t;
+		  sinix *= exp_t;
+		  cosix *= exp_t;
+		}
+	    }
+	  if (__real__ x > t)
+	    {
+	      /* Overflow (original real part of x > 3t).  */
+	      __real__ retval = DBL_MAX * cosix;
+	      __imag__ retval = DBL_MAX * sinix;
 	    }
 	  else
 	    {
-	      __real__ retval = __copysign (exp_val, cosix);
-	      __imag__ retval = __copysign (exp_val, sinix);
+	      double exp_val = __ieee754_exp (__real__ x);
+	      __real__ retval = exp_val * cosix;
+	      __imag__ retval = exp_val * sinix;
 	    }
 	}
       else
diff --git a/math/s_cexpf.c b/math/s_cexpf.c
index a9c28ed..4aa9765 100644
--- a/math/s_cexpf.c
+++ b/math/s_cexpf.c
@@ -1,5 +1,5 @@
 /* Return value of complex exponential function for float complex value.
-   Copyright (C) 1997, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -21,7 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ float
 __cexpf (__complex__ float x)
@@ -36,20 +36,35 @@ __cexpf (__complex__ float x)
       if (__builtin_expect (icls >= FP_ZERO, 1))
 	{
 	  /* Imaginary part is finite.  */
-	  float exp_val = __ieee754_expf (__real__ x);
+	  const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2);
 	  float sinix, cosix;
 
 	  __sincosf (__imag__ x, &sinix, &cosix);
 
-	  if (isfinite (exp_val))
+	  if (__real__ x > t)
 	    {
-	      __real__ retval = exp_val * cosix;
-	      __imag__ retval = exp_val * sinix;
+	      float exp_t = __ieee754_expf (t);
+	      __real__ x -= t;
+	      sinix *= exp_t;
+	      cosix *= exp_t;
+	      if (__real__ x > t)
+		{
+		  __real__ x -= t;
+		  sinix *= exp_t;
+		  cosix *= exp_t;
+		}
+	    }
+	  if (__real__ x > t)
+	    {
+	      /* Overflow (original real part of x > 3t).  */
+	      __real__ retval = FLT_MAX * cosix;
+	      __imag__ retval = FLT_MAX * sinix;
 	    }
 	  else
 	    {
-	      __real__ retval = __copysignf (exp_val, cosix);
-	      __imag__ retval = __copysignf (exp_val, sinix);
+	      float exp_val = __ieee754_expf (__real__ x);
+	      __real__ retval = exp_val * cosix;
+	      __imag__ retval = exp_val * sinix;
 	    }
 	}
       else
diff --git a/math/s_cexpl.c b/math/s_cexpl.c
index 3059880..2568249 100644
--- a/math/s_cexpl.c
+++ b/math/s_cexpl.c
@@ -1,5 +1,5 @@
 /* Return value of complex exponential function for long double complex value.
-   Copyright (C) 1997, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -21,7 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ long double
 __cexpl (__complex__ long double x)
@@ -36,20 +36,35 @@ __cexpl (__complex__ long double x)
       if (__builtin_expect (icls >= FP_ZERO, 1))
 	{
 	  /* Imaginary part is finite.  */
-	  long double exp_val = __ieee754_expl (__real__ x);
+	  const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l);
 	  long double sinix, cosix;
 
 	  __sincosl (__imag__ x, &sinix, &cosix);
 
-	  if (isfinite (exp_val))
+	  if (__real__ x > t)
 	    {
-	      __real__ retval = exp_val * cosix;
-	      __imag__ retval = exp_val * sinix;
+	      long double exp_t = __ieee754_expl (t);
+	      __real__ x -= t;
+	      sinix *= exp_t;
+	      cosix *= exp_t;
+	      if (__real__ x > t)
+		{
+		  __real__ x -= t;
+		  sinix *= exp_t;
+		  cosix *= exp_t;
+		}
+	    }
+	  if (__real__ x > t)
+	    {
+	      /* Overflow (original real part of x > 3t).  */
+	      __real__ retval = LDBL_MAX * cosix;
+	      __imag__ retval = LDBL_MAX * sinix;
 	    }
 	  else
 	    {
-	      __real__ retval = __copysignl (exp_val, cosix);
-	      __imag__ retval = __copysignl (exp_val, sinix);
+	      long double exp_val = __ieee754_expl (__real__ x);
+	      __real__ retval = exp_val * cosix;
+	      __imag__ retval = exp_val * sinix;
 	    }
 	}
       else