NPTL is no longer an add-on!

[platform/upstream/glibc.git] / math / s_ctan.c

diff --git a/math/s_ctan.c b/math/s_ctan.c
index c838fad..407ff5a 100644 (file)
--- a/math/s_ctan.c
+++ b/math/s_ctan.c
@@ -1,5 +1,5 @@
 /* Complex tangent function for double.
-   Copyright (C) 1997, 2005, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2014 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -20,16 +20,15 @@
 #include <complex.h>
 #include <fenv.h>
 #include <math.h>
-
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ double
 __ctan (__complex__ double x)
 {
   __complex__ double res;
 
-  if (__builtin_expect (!isfinite (__real__ x) || !isfinite (__imag__ x), 0))
+  if (__glibc_unlikely (!isfinite (__real__ x) || !isfinite (__imag__ x)))
     {
       if (__isinf_ns (__imag__ x))
        {
@@ -51,24 +50,66 @@ __ctan (__complex__ double x)
     }
   else
     {
-      double sin2rx, cos2rx;
+      double sinrx, cosrx;
       double den;
+      const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2 / 2);
+      int rcls = fpclassify (__real__ x);
 
-      __sincos (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_cosh (2.0 * __imag__ x);
+      if (__glibc_likely (rcls != FP_SUBNORMAL))
+       {
+         __sincos (__real__ x, &sinrx, &cosrx);
+       }
+      else
+       {
+         sinrx = __real__ x;
+         cosrx = 1.0;
+       }
 
-      if (den == 0.0)
+      if (fabs (__imag__ x) > t)
        {
-         __complex__ double ez = __cexp (1.0i * x);
-         __complex__ double emz = __cexp (-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).  */
+         double exp_2t = __ieee754_exp (2 * t);
 
-         res = (ez - emz) / (ez + emz) * -1.0i;
+         __imag__ res = __copysign (1.0, __imag__ x);
+         __real__ res = 4 * sinrx * cosrx;
+         __imag__ x = fabs (__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_exp (2 * __imag__ x);
        }
       else
        {
-         __real__ res = sin2rx / den;
-         __imag__ res = __ieee754_sinh (2.0 * __imag__ x) / den;
+         double sinhix, coshix;
+         if (fabs (__imag__ x) > DBL_MIN)
+           {
+             sinhix = __ieee754_sinh (__imag__ x);
+             coshix = __ieee754_cosh (__imag__ x);
+           }
+         else
+           {
+             sinhix = __imag__ x;
+             coshix = 1.0;
+           }
+
+         if (fabs (sinhix) > fabs (cosrx) * DBL_EPSILON)
+           den = cosrx * cosrx + sinhix * sinhix;
+         else
+           den = cosrx * cosrx;
+         __real__ res = sinrx * cosrx / den;
+         __imag__ res = sinhix * coshix / den;
        }
     }