diff options
Diffstat (limited to 'math/k_casinh.c')
| -rw-r--r-- | math/k_casinh.c | 32 |
1 files changed, 32 insertions, 0 deletions
diff --git a/math/k_casinh.c b/math/k_casinh.c index 45845e2164..4cb232a17a 100644 --- a/math/k_casinh.c +++ b/math/k_casinh.c @@ -77,6 +77,38 @@ __kernel_casinh (__complex__ double x, int adj) else __imag__ res = __ieee754_atan2 (s, rx); } + else if (ix > 1.0 && ix < 1.5 && rx < 0.5) + { + if (rx < DBL_EPSILON * DBL_EPSILON) + { + double ix2m1 = (ix + 1.0) * (ix - 1.0); + double s = __ieee754_sqrt (ix2m1); + + __real__ res = __log1p (2.0 * (ix2m1 + ix * s)) / 2.0; + if (adj) + __imag__ res = __ieee754_atan2 (rx, __copysign (s, __imag__ x)); + else + __imag__ res = __ieee754_atan2 (s, rx); + } + else + { + double ix2m1 = (ix + 1.0) * (ix - 1.0); + double rx2 = rx * rx; + double f = rx2 * (2.0 + rx2 + 2.0 * ix * ix); + double d = __ieee754_sqrt (ix2m1 * ix2m1 + f); + double dp = d + ix2m1; + double dm = f / dp; + double r1 = __ieee754_sqrt ((dm + rx2) / 2.0); + double r2 = rx * ix / r1; + + __real__ res = __log1p (rx2 + dp + 2.0 * (rx * r1 + ix * r2)) / 2.0; + if (adj) + __imag__ res = __ieee754_atan2 (rx + r1, __copysign (ix + r2, + __imag__ x)); + else + __imag__ res = __ieee754_atan2 (ix + r2, rx + r1); + } + } else if (ix == 1.0 && rx < 0.5) { if (rx < DBL_EPSILON / 8.0) |