Fix cexp overflow (bug 13892).

1 files changed, 23 insertions, 8 deletions

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