Fix ctan, ctanh overflow (bug 11521).

1 files changed, 32 insertions, 11 deletions

diff --git a/math/s_ctanhf.c b/math/s_ctanhf.c
index fce5aaf..e505155 100644
--- a/math/s_ctanhf.c
+++ b/math/s_ctanhf.c
@@ -1,5 +1,5 @@
 /* Complex hyperbole tangent for float.
-   Copyright (C) 1997, 2005, 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
 __ctanhf (__complex__ float x)
@@ -50,24 +50,45 @@ __ctanhf (__complex__ float x)
     }
   else
     {
-      float sin2ix, cos2ix;
+      float sinix, cosix;
       float den;
+      const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2 / 2);
 
-      __sincosf (2.0 * __imag__ x, &sin2ix, &cos2ix);
+      /* tanh(x+iy) = (sinh(2x) + i*sin(2y))/(cosh(2x) + cos(2y))
+	 = (sinh(x)*cosh(x) + i*sin(y)*cos(y))/(sinh(x)^2 + cos(y)^2).  */
 
-      den = (__ieee754_coshf (2.0 * __real__ x) + cos2ix);
+      __sincosf (__imag__ x, &sinix, &cosix);
 
-      if (den == 0.0f)
+      if (fabsf (__real__ x) > t)
 	{
-	  __complex__ float ez = __cexpf (x);
-	  __complex__ float emz = __cexpf (-x);
+	  /* Avoid intermediate overflow when the imaginary part of
+	     the result may be subnormal.  Ignoring negligible terms,
+	     the real part is +/- 1, the imaginary part is
+	     sin(y)*cos(y)/sinh(x)^2 = 4*sin(y)*cos(y)/exp(2x).  */
+	  float exp_2t = __ieee754_expf (2 * t);
 
-	  res = (ez - emz) / (ez + emz);
+	  __real__ res = __copysignf (1.0, __real__ x);
+	  __imag__ res = 4 * sinix * cosix;
+	  __real__ x = fabsf (__real__ x);
+	  __real__ x -= t;
+	  __imag__ res /= exp_2t;
+	  if (__real__ x > t)
+	    {
+	      /* Underflow (original real part of x has absolute value
+		 > 2t).  */
+	      __imag__ res /= exp_2t;
+	    }
+	  else
+	    __imag__ res /= __ieee754_expf (2 * __real__ x);
 	}
       else
 	{
-	  __real__ res = __ieee754_sinhf (2.0 * __real__ x) / den;
-	  __imag__ res = sin2ix / den;
+	  float sinhrx = __ieee754_sinhf (__real__ x);
+	  float coshrx = __ieee754_coshf (__real__ x);
+
+	  den = sinhrx * sinhrx + cosix * cosix;
+	  __real__ res = sinhrx * coshrx / den;
+	  __imag__ res = sinix * cosix / den;
 	}
     }