diff options
Diffstat (limited to 'math/s_ctanl.c')
| -rw-r--r-- | math/s_ctanl.c | 45 |
1 files changed, 32 insertions, 13 deletions
diff --git a/math/s_ctanl.c b/math/s_ctanl.c index 112dd723d8..12d700cad9 100644 --- a/math/s_ctanl.c +++ b/math/s_ctanl.c @@ -1,5 +1,5 @@ /* Complex tangent function for long double. - 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. @@ -20,9 +20,8 @@ #include <complex.h> #include <fenv.h> #include <math.h> - #include <math_private.h> - +#include <float.h> __complex__ long double __ctanl (__complex__ long double x) @@ -51,25 +50,45 @@ __ctanl (__complex__ long double x) } 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); - - den = cos2rx + __ieee754_coshl (2.0 * __imag__ x); + /* 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). */ + __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; } } |