From bcc8d6617ba029c288fff9680a02b9a3b1caa9c0 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Mon, 9 Apr 2012 22:31:35 +0000 Subject: Fix ctan, ctanh overflow (bug 11521). --- math/s_ctanhl.c | 43 ++++++++++++++++++++++++++++++++----------- 1 file changed, 32 insertions(+), 11 deletions(-) (limited to 'math/s_ctanhl.c') diff --git a/math/s_ctanhl.c b/math/s_ctanhl.c index d22e13a..e5d6779 100644 --- a/math/s_ctanhl.c +++ b/math/s_ctanhl.c @@ -1,5 +1,5 @@ /* Complex hyperbole tangent 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 , 1997. @@ -21,7 +21,7 @@ #include #include #include - +#include __complex__ long double __ctanhl (__complex__ long double x) @@ -50,24 +50,45 @@ __ctanhl (__complex__ long double x) } else { - long double sin2ix, cos2ix; + long double sinix, cosix; long double den; + const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2); - __sincosl (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_coshl (2.0 * __real__ x) + cos2ix); + __sincosl (__imag__ x, &sinix, &cosix); - if (den == 0.0L) + if (fabsl (__real__ x) > t) { - __complex__ long double ez = __cexpl (x); - __complex__ long double emz = __cexpl (-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). */ + long double exp_2t = __ieee754_expl (2 * t); - res = (ez - emz) / (ez + emz); + __real__ res = __copysignl (1.0, __real__ x); + __imag__ res = 4 * sinix * cosix; + __real__ x = fabsl (__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_expl (2 * __real__ x); } else { - __real__ res = __ieee754_sinhl (2.0 * __real__ x) / den; - __imag__ res = sin2ix / den; + long double sinhrx = __ieee754_sinhl (__real__ x); + long double coshrx = __ieee754_coshl (__real__ x); + + den = sinhrx * sinhrx + cosix * cosix; + __real__ res = sinhrx * coshrx / den; + __imag__ res = sinix * cosix / den; } } -- cgit v1.1