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_ctanhf.c | 43 ++++++++++++++++++++++++++++++++----------- 1 file changed, 32 insertions(+), 11 deletions(-) (limited to 'math/s_ctanhf.c') 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 , 1997. @@ -21,7 +21,7 @@ #include #include #include - +#include __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; } } -- cgit v1.1