aboutsummaryrefslogtreecommitdiff
path: root/math/s_ctanhl.c
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2012-04-09 22:31:35 +0000
committerJoseph Myers <joseph@codesourcery.com>2012-04-09 22:31:35 +0000
commitbcc8d6617ba029c288fff9680a02b9a3b1caa9c0 (patch)
tree9cb163731d4d165bfe71bad55f59d88c6b420917 /math/s_ctanhl.c
parent8a1fbaaf75536088b4d505409e0e87975d8cf8d5 (diff)
downloadglibc-bcc8d6617ba029c288fff9680a02b9a3b1caa9c0.zip
glibc-bcc8d6617ba029c288fff9680a02b9a3b1caa9c0.tar.gz
glibc-bcc8d6617ba029c288fff9680a02b9a3b1caa9c0.tar.bz2
Fix ctan, ctanh overflow (bug 11521).
Diffstat (limited to 'math/s_ctanhl.c')
-rw-r--r--math/s_ctanhl.c43
1 files changed, 32 insertions, 11 deletions
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 <drepper@cygnus.com>, 1997.
@@ -21,7 +21,7 @@
#include <fenv.h>
#include <math.h>
#include <math_private.h>
-
+#include <float.h>
__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;
}
}