diff options
Diffstat (limited to 'math/s_catanhl.c')
| -rw-r--r-- | math/s_catanhl.c | 86 |
1 files changed, 78 insertions, 8 deletions
diff --git a/math/s_catanhl.c b/math/s_catanhl.c index 6844f03dea..8c4b8940c0 100644 --- a/math/s_catanhl.c +++ b/math/s_catanhl.c @@ -20,7 +20,14 @@ #include <complex.h> #include <math.h> #include <math_private.h> +#include <float.h> +/* To avoid spurious overflows, use this definition to treat IBM long + double as approximating an IEEE-style format. */ +#if LDBL_MANT_DIG == 106 +# undef LDBL_EPSILON +# define LDBL_EPSILON 0x1p-106L +#endif __complex__ long double __catanhl (__complex__ long double x) @@ -56,19 +63,82 @@ __catanhl (__complex__ long double x) } else { - long double i2 = __imag__ x * __imag__ x; + if (fabsl (__real__ x) >= 16.0L / LDBL_EPSILON + || fabsl (__imag__ x) >= 16.0L / LDBL_EPSILON) + { + __imag__ res = __copysignl (M_PI_2l, __imag__ x); + if (fabsl (__imag__ x) <= 1.0L) + __real__ res = 1.0L / __real__ x; + else if (fabsl (__real__ x) <= 1.0L) + __real__ res = __real__ x / __imag__ x / __imag__ x; + else + { + long double h = __ieee754_hypotl (__real__ x / 2.0L, + __imag__ x / 2.0L); + __real__ res = __real__ x / h / h / 4.0L; + } + } + else + { + if (fabsl (__real__ x) == 1.0L + && fabsl (__imag__ x) < LDBL_EPSILON * LDBL_EPSILON) + __real__ res = (__copysignl (0.5L, __real__ x) + * (M_LN2l - __ieee754_logl (fabsl (__imag__ x)))); + else + { + long double i2 = 0.0; + if (fabsl (__imag__ x) >= LDBL_EPSILON * LDBL_EPSILON) + i2 = __imag__ x * __imag__ x; + + long double num = 1.0L + __real__ x; + num = i2 + num * num; - long double num = 1.0 + __real__ x; - num = i2 + num * num; + long double den = 1.0L - __real__ x; + den = i2 + den * den; - long double den = 1.0 - __real__ x; - den = i2 + den * den; + long double f = num / den; + if (f < 0.5L) + __real__ res = 0.25L * __ieee754_logl (f); + else + { + num = 4.0L * __real__ x; + __real__ res = 0.25L * __log1pl (num / den); + } + } - __real__ res = 0.25 * (__ieee754_logl (num) - __ieee754_logl (den)); + long double absx, absy, den; - den = 1 - __real__ x * __real__ x - i2; + absx = fabsl (__real__ x); + absy = fabsl (__imag__ x); + if (absx < absy) + { + long double t = absx; + absx = absy; + absy = t; + } - __imag__ res = 0.5 * __ieee754_atan2l (2.0 * __imag__ x, den); + if (absy < LDBL_EPSILON / 2.0L) + den = (1.0L - absx) * (1.0L + absx); + else if (absx >= 1.0L) + den = (1.0L - absx) * (1.0L + absx) - absy * absy; + else if (absx >= 0.75L || absy >= 0.5L) + den = -__x2y2m1l (absx, absy); + else + den = (1.0L - absx) * (1.0L + absx) - absy * absy; + + __imag__ res = 0.5L * __ieee754_atan2l (2.0L * __imag__ x, den); + } + + if (fabsl (__real__ res) < LDBL_MIN) + { + volatile long double force_underflow = __real__ res * __real__ res; + (void) force_underflow; + } + if (fabsl (__imag__ res) < LDBL_MIN) + { + volatile long double force_underflow = __imag__ res * __imag__ res; + (void) force_underflow; + } } return res; |