Fix ctan, ctanh overflow for ldbl-128ibm (bug 11521).

1 files changed, 34 insertions, 14 deletions

diff --git a/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c b/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c
index 0b1bc4f70e..9d89bbe311 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c
@@ -1,5 +1,5 @@
 /* Complex tangent function for long double.  IBM extended format version.
-   Copyright (C) 1997,2005,2006 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,6 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_ldbl_opt.h>
+#include <float.h>
 
 #include <math_private.h>
 
@@ -46,34 +47,53 @@ __ctanl (__complex__ long double x)
 	  __real__ res = __nanl ("");
 	  __imag__ res = __nanl ("");
 
-#ifdef FE_INVALID
-	  if (__isinfl (__real__ x))
+	  if (__isinf_nsl (__real__ x))
 	    feraiseexcept (FE_INVALID);
-#endif
 	}
     }
   else
     {
-      long double sin2rx, cos2rx;
+      long double sinrx, cosrx;
       long double den;
+      const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2);
 
-      __sincosl (2.0 * __real__ x, &sin2rx, &cos2rx);
+      /* tan(x+iy) = (sin(2x) + i*sinh(2y))/(cos(2x) + cosh(2y))
+        = (sin(x)*cos(x) + i*sinh(y)*cosh(y)/(cos(x)^2 + sinh(y)^2). */
 
-      den = cos2rx + __ieee754_coshl (2.0 * __imag__ x);
+      __sincosl (__real__ x, &sinrx, &cosrx);
 
-
-      if (den == 0.0)
+      if (fabsl (__imag__ x) > t)
 	{
-	  __complex__ long double ez = __cexpl (1.0i * x);
-	  __complex__ long double emz = __cexpl (-1.0i * x);
+	  /* Avoid intermediate overflow when the real part of the
+	    result may be subnormal.  Ignoring negligible terms, the
+	    imaginary part is +/- 1, the real part is
+	    sin(x)*cos(x)/sinh(y)^2 = 4*sin(x)*cos(x)/exp(2y).  */
+	  long double exp_2t = __ieee754_expl (2 * t);
 
-	  res = (ez - emz) / (ez + emz) * -1.0i;
+	  __imag__ res = __copysignl (1.0, __imag__ x);
+	  __real__ res = 4 * sinrx * cosrx;
+	  __imag__ x = fabsl (__imag__ x);
+	  __imag__ x -= t;
+	  __real__ res /= exp_2t;
+	  if (__imag__ x > t)
+	    {
+	      /* Underflow (original imaginary part of x has absolute
+		 value > 2t).  */
+	      __real__ res /= exp_2t;
+	    }
+	  else
+	    __real__ res /= __ieee754_expl (2 * __imag__ x);
 	}
       else
 	{
-	  __real__ res = sin2rx / den;
-	  __imag__ res = __ieee754_sinhl (2.0 * __imag__ x) / den;
+	  long double sinhix = __ieee754_sinhl (__imag__ x);
+	  long double coshix = __ieee754_coshl (__imag__ x);
+
+	  den = cosrx * cosrx + sinhix * sinhix;
+	  __real__ res = sinrx * cosrx / den;
+	  __imag__ res = sinhix * coshix / den;
 	}
+
       /* __gcc_qmul does not respect -0.0 so we need the following fixup.  */
       if ((__real__ res == 0.0) && (__real__ x == 0.0))
         __real__ res = __real__ x;