Fix casinh inaccuracy for imaginary part < 1.0, real part small (bug 10357).

1 files changed, 52 insertions, 0 deletions

diff --git a/math/k_casinhl.c b/math/k_casinhl.c
index 8e5bbd4385..b5f94006f7 100644
--- a/math/k_casinhl.c
+++ b/math/k_casinhl.c
@@ -143,6 +143,58 @@ __kernel_casinhl (__complex__ long double x, int adj)
 	    __imag__ res = __ieee754_atan2l (1.0L + s2, rx + s1);
 	}
     }
+  else if (ix < 1.0L && rx < 0.5L)
+    {
+      if (ix >= LDBL_EPSILON)
+	{
+	  if (rx < LDBL_EPSILON * LDBL_EPSILON)
+	    {
+	      long double onemix2 = (1.0L + ix) * (1.0L - ix);
+	      long double s = __ieee754_sqrtl (onemix2);
+
+	      __real__ res = __log1pl (2.0L * rx / s) / 2.0L;
+	      if (adj)
+		__imag__ res = __ieee754_atan2l (s, __imag__ x);
+	      else
+		__imag__ res = __ieee754_atan2l (ix, s);
+	    }
+	  else
+	    {
+	      long double onemix2 = (1.0L + ix) * (1.0L - ix);
+	      long double rx2 = rx * rx;
+	      long double f = rx2 * (2.0L + rx2 + 2.0L * ix * ix);
+	      long double d = __ieee754_sqrtl (onemix2 * onemix2 + f);
+	      long double dp = d + onemix2;
+	      long double dm = f / dp;
+	      long double r1 = __ieee754_sqrtl ((dp + rx2) / 2.0L);
+	      long double r2 = rx * ix / r1;
+
+	      __real__ res
+		= __log1pl (rx2 + dm + 2.0L * (rx * r1 + ix * r2)) / 2.0L;
+	      if (adj)
+		__imag__ res = __ieee754_atan2l (rx + r1,
+						 __copysignl (ix + r2,
+							      __imag__ x));
+	      else
+		__imag__ res = __ieee754_atan2l (ix + r2, rx + r1);
+	    }
+	}
+      else
+	{
+	  long double s = __ieee754_hypotl (1.0L, rx);
+
+	  __real__ res = __log1pl (2.0L * rx * (rx + s)) / 2.0L;
+	  if (adj)
+	    __imag__ res = __ieee754_atan2l (s, __imag__ x);
+	  else
+	    __imag__ res = __ieee754_atan2l (ix, s);
+	}
+      if (__real__ res < LDBL_MIN)
+	{
+	  volatile long double force_underflow = __real__ res * __real__ res;
+	  (void) force_underflow;
+	}
+    }
   else
     {
       __real__ y = (rx - ix) * (rx + ix) + 1.0L;