Fix ccos, csin, ccosh, csinh overflows (bug 14123).

1 files changed, 36 insertions, 7 deletions

diff --git a/math/s_ccoshf.c b/math/s_ccoshf.c
index a5a12c1..d6f8110 100644
--- a/math/s_ccoshf.c
+++ b/math/s_ccoshf.c
@@ -1,5 +1,5 @@
 /* Complex cosine hyperbole function for float.
-   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.
 
@@ -20,9 +20,8 @@
 #include <complex.h>
 #include <fenv.h>
 #include <math.h>
-
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ float
 __ccoshf (__complex__ float x)
@@ -37,14 +36,44 @@ __ccoshf (__complex__ float x)
       if (__builtin_expect (icls >= FP_ZERO, 1))
 	{
 	  /* Imaginary part is finite.  */
-	  float sinh_val = __ieee754_sinhf (__real__ x);
-	  float cosh_val = __ieee754_coshf (__real__ x);
+	  const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2);
 	  float sinix, cosix;
 
 	  __sincosf (__imag__ x, &sinix, &cosix);
 
-	  __real__ retval = cosh_val * cosix;
-	  __imag__ retval = sinh_val * sinix;
+	  if (fabsf (__real__ x) > t)
+	    {
+	      float exp_t = __ieee754_expf (t);
+	      float rx = fabsf (__real__ x);
+	      if (signbit (__real__ x))
+		sinix = -sinix;
+	      rx -= t;
+	      sinix *= exp_t / 2.0f;
+	      cosix *= exp_t / 2.0f;
+	      if (rx > t)
+		{
+		  rx -= t;
+		  sinix *= exp_t;
+		  cosix *= exp_t;
+		}
+	      if (rx > t)
+		{
+		  /* Overflow (original real part of x > 3t).  */
+		  __real__ retval = FLT_MAX * cosix;
+		  __imag__ retval = FLT_MAX * sinix;
+		}
+	      else
+		{
+		  float exp_val = __ieee754_expf (rx);
+		  __real__ retval = exp_val * cosix;
+		  __imag__ retval = exp_val * sinix;
+		}
+	    }
+	  else
+	    {
+	      __real__ retval = __ieee754_coshf (__real__ x) * cosix;
+	      __imag__ retval = __ieee754_sinhf (__real__ x) * sinix;
+	    }
 	}
       else
 	{