Fix cacos real-part inaccuracy for result real part near 0 (bug 15023).

16 files changed, 430 insertions, 130 deletions

diff --git a/ChangeLog b/ChangeLog
index e7cfb64231..890e3d4930 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,34 @@
+2013-01-17  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #15023]
+	* include/complex.h: Condition contents on [!_COMPLEX_H].
+	(__kernel_casinhf): New prototype.
+	(__kernel_casinh): Likewise.
+	(__kernel_casinhl): Likewise.
+	* math/Makefile (libm_calls): Add k_casinh.
+	* math/k_casinh.c: New file.
+	* math/k_casinhf.c: Likewise.
+	* math/k_casinhl.c: Likewise.
+	* math/s_cacos.c (__cacos): Implement using __kernel_casinh for
+	finite nonzero arguments.
+	* math/s_cacosf.c (__cacosf): Implement using __kernel_casinhf for
+	finite nonzero arguments.
+	* math/s_cacosl.c (__cacosl): Implement using __kernel_casinhl for
+	finite nonzero arguments.
+	* math/s_casinh.c: Do not include <float.h>.
+	(__casinh): Move code for finite nonzero arguments to k_casinh.c.
+	* math/s_casinhf.c: Do not include <float.h>.
+	(__casinhf): Move code for finite nonzero arguments to
+	k_casinhf.c.
+	* math/s_casinhl.c: Do not include <float.h>.
+	[LDBL_MANT_DIG == 106] (LDBL_EPSILON): Do not undefine and
+	redefine.
+	(__casinhl): Move code for finite nonzero arguments to
+	k_casinhl.c.
+	* math/libm-test.inc (cacos_test): Add more tests.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+
 2013-01-17  Pino Toscano  <toscano.pino@tiscali.it>
 
 	* sysdeps/unix/sysv/linux/malloc-sysdep.h (HAVE_MREMAP): New define.
diff --git a/NEWS b/NEWS
index 55f65fd686..99fb54226f 100644
--- a/NEWS
+++ b/NEWS
@@ -10,7 +10,7 @@ Version 2.18
 * The following bugs are resolved with this release:
 
   13951, 14200, 14317, 14327, 14964, 14981, 14982, 14985, 14994, 14996,
-  15003.
+  15003, 15023.
 
 
 Version 2.17
diff --git a/include/complex.h b/include/complex.h
index acf8cf14ba..e173f1f6a3 100644
--- a/include/complex.h
+++ b/include/complex.h
@@ -1 +1,11 @@
-#include <math/complex.h>
+#ifndef _COMPLEX_H
+# include <math/complex.h>
+
+/* Return the complex inverse hyperbolic sine of finite nonzero Z,
+   with the imaginary part of the result subtracted from pi/2 if ADJ
+   is nonzero.  */
+extern complex float __kernel_casinhf (complex float z, int adj);
+extern complex double __kernel_casinh (complex double z, int adj);
+extern complex long double __kernel_casinhl (complex long double z, int adj);
+
+#endif
diff --git a/math/Makefile b/math/Makefile
index b9519cfc24..da18b56d4b 100644
--- a/math/Makefile
+++ b/math/Makefile
@@ -58,7 +58,7 @@ libm-calls = e_acos e_acosh e_asin e_atan2 e_atanh e_cosh e_exp e_fmod	\
 	     s_catan s_casin s_ccos s_csin s_ctan s_ctanh s_cacos	\
 	     s_casinh s_cacosh s_catanh s_csqrt s_cpow s_cproj s_clog10 \
 	     s_fma s_lrint s_llrint s_lround s_llround e_exp10 w_log2	\
-	     s_isinf_ns $(calls:s_%=m_%) x2y2m1
+	     s_isinf_ns $(calls:s_%=m_%) x2y2m1 k_casinh
 
 include ../Makeconfig
 
diff --git a/math/k_casinh.c b/math/k_casinh.c
new file mode 100644
index 0000000000..7f98f24a80
--- /dev/null
+++ b/math/k_casinh.c
@@ -0,0 +1,85 @@
+/* Return arc hyperbole sine for double value, with the imaginary part
+   of the result possibly adjusted for use in computing other
+   functions.
+   Copyright (C) 1997-2013 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <complex.h>
+#include <math.h>
+#include <math_private.h>
+#include <float.h>
+
+/* Return the complex inverse hyperbolic sine of finite nonzero Z,
+   with the imaginary part of the result subtracted from pi/2 if ADJ
+   is nonzero.  */
+
+__complex__ double
+__kernel_casinh (__complex__ double x, int adj)
+{
+  __complex__ double res;
+  double rx, ix;
+  __complex__ double y;
+
+  /* Avoid cancellation by reducing to the first quadrant.  */
+  rx = fabs (__real__ x);
+  ix = fabs (__imag__ x);
+
+  if (rx >= 1.0 / DBL_EPSILON || ix >= 1.0 / DBL_EPSILON)
+    {
+      /* For large x in the first quadrant, x + csqrt (1 + x * x)
+	 is sufficiently close to 2 * x to make no significant
+	 difference to the result; avoid possible overflow from
+	 the squaring and addition.  */
+      __real__ y = rx;
+      __imag__ y = ix;
+
+      if (adj)
+	{
+	  double t = __real__ y;
+	  __real__ y = __copysign (__imag__ y, __imag__ x);
+	  __imag__ y = t;
+	}
+
+      res = __clog (y);
+      __real__ res += M_LN2;
+    }
+  else
+    {
+      __real__ y = (rx - ix) * (rx + ix) + 1.0;
+      __imag__ y = 2.0 * rx * ix;
+
+      y = __csqrt (y);
+
+      __real__ y += rx;
+      __imag__ y += ix;
+
+      if (adj)
+	{
+	  double t = __real__ y;
+	  __real__ y = copysign (__imag__ y, __imag__ x);
+	  __imag__ y = t;
+	}
+
+      res = __clog (y);
+    }
+
+  /* Give results the correct sign for the original argument.  */
+  __real__ res = __copysign (__real__ res, __real__ x);
+  __imag__ res = __copysign (__imag__ res, (adj ? 1.0 : __imag__ x));
+
+  return res;
+}
diff --git a/math/k_casinhf.c b/math/k_casinhf.c
new file mode 100644
index 0000000000..9401636348
--- /dev/null
+++ b/math/k_casinhf.c
@@ -0,0 +1,85 @@
+/* Return arc hyperbole sine for float value, with the imaginary part
+   of the result possibly adjusted for use in computing other
+   functions.
+   Copyright (C) 1997-2013 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <complex.h>
+#include <math.h>
+#include <math_private.h>
+#include <float.h>
+
+/* Return the complex inverse hyperbolic sine of finite nonzero Z,
+   with the imaginary part of the result subtracted from pi/2 if ADJ
+   is nonzero.  */
+
+__complex__ float
+__kernel_casinhf (__complex__ float x, int adj)
+{
+  __complex__ float res;
+  float rx, ix;
+  __complex__ float y;
+
+  /* Avoid cancellation by reducing to the first quadrant.  */
+  rx = fabsf (__real__ x);
+  ix = fabsf (__imag__ x);
+
+  if (rx >= 1.0f / FLT_EPSILON || ix >= 1.0f / FLT_EPSILON)
+    {
+      /* For large x in the first quadrant, x + csqrt (1 + x * x)
+	 is sufficiently close to 2 * x to make no significant
+	 difference to the result; avoid possible overflow from
+	 the squaring and addition.  */
+      __real__ y = rx;
+      __imag__ y = ix;
+
+      if (adj)
+	{
+	  float t = __real__ y;
+	  __real__ y = __copysignf (__imag__ y, __imag__ x);
+	  __imag__ y = t;
+	}
+
+      res = __clogf (y);
+      __real__ res += (float) M_LN2;
+    }
+  else
+    {
+      __real__ y = (rx - ix) * (rx + ix) + 1.0;
+      __imag__ y = 2.0 * rx * ix;
+
+      y = __csqrtf (y);
+
+      __real__ y += rx;
+      __imag__ y += ix;
+
+      if (adj)
+	{
+	  float t = __real__ y;
+	  __real__ y = __copysignf (__imag__ y, __imag__ x);
+	  __imag__ y = t;
+	}
+
+      res = __clogf (y);
+    }
+
+  /* Give results the correct sign for the original argument.  */
+  __real__ res = __copysignf (__real__ res, __real__ x);
+  __imag__ res = __copysignf (__imag__ res, (adj ? 1.0f : __imag__ x));
+
+  return res;
+}
diff --git a/math/k_casinhl.c b/math/k_casinhl.c
new file mode 100644
index 0000000000..6412979755
--- /dev/null
+++ b/math/k_casinhl.c
@@ -0,0 +1,92 @@
+/* Return arc hyperbole sine for long double value, with the imaginary
+   part of the result possibly adjusted for use in computing other
+   functions.
+   Copyright (C) 1997-2013 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#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
+
+/* Return the complex inverse hyperbolic sine of finite nonzero Z,
+   with the imaginary part of the result subtracted from pi/2 if ADJ
+   is nonzero.  */
+
+__complex__ long double
+__kernel_casinhl (__complex__ long double x, int adj)
+{
+  __complex__ long double res;
+  long double rx, ix;
+  __complex__ long double y;
+
+  /* Avoid cancellation by reducing to the first quadrant.  */
+  rx = fabsl (__real__ x);
+  ix = fabsl (__imag__ x);
+
+  if (rx >= 1.0L / LDBL_EPSILON || ix >= 1.0L / LDBL_EPSILON)
+    {
+      /* For large x in the first quadrant, x + csqrt (1 + x * x)
+	 is sufficiently close to 2 * x to make no significant
+	 difference to the result; avoid possible overflow from
+	 the squaring and addition.  */
+      __real__ y = rx;
+      __imag__ y = ix;
+
+      if (adj)
+	{
+	  long double t = __real__ y;
+	  __real__ y = __copysignl (__imag__ y, __imag__ x);
+	  __imag__ y = t;
+	}
+
+      res = __clogl (y);
+      __real__ res += M_LN2l;
+    }
+  else
+    {
+      __real__ y = (rx - ix) * (rx + ix) + 1.0;
+      __imag__ y = 2.0 * rx * ix;
+
+      y = __csqrtl (y);
+
+      __real__ y += rx;
+      __imag__ y += ix;
+
+      if (adj)
+	{
+	  long double t = __real__ y;
+	  __real__ y = __copysignl (__imag__ y, __imag__ x);
+	  __imag__ y = t;
+	}
+
+      res = __clogl (y);
+    }
+
+  /* Give results the correct sign for the original argument.  */
+  __real__ res = __copysignl (__real__ res, __real__ x);
+  __imag__ res = __copysignl (__imag__ res, (adj ? 1.0L : __imag__ x));
+
+  return res;
+}
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 56e321781e..1c2970fc79 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -1453,6 +1453,43 @@ cacos_test (void)
   TEST_c_c (cacos, 1.5L, plus_zero, plus_zero, -0.9624236501192068949955178268487368462704L);
   TEST_c_c (cacos, 1.5L, minus_zero, plus_zero, 0.9624236501192068949955178268487368462704L);
 
+  TEST_c_c (cacos, 0x1p50L, 1.0L, 8.881784197001252323389053344727730248720e-16L, -3.535050620855721078027883819436720218708e1L);
+  TEST_c_c (cacos, 0x1p50L, -1.0L, 8.881784197001252323389053344727730248720e-16L, 3.535050620855721078027883819436720218708e1L);
+  TEST_c_c (cacos, -0x1p50L, 1.0L, 3.141592653589792350284223683154270545292L, -3.535050620855721078027883819436720218708e1L);
+  TEST_c_c (cacos, -0x1p50L, -1.0L, 3.141592653589792350284223683154270545292L, 3.535050620855721078027883819436720218708e1L);
+  TEST_c_c (cacos, 1.0L, 0x1p50L, 1.570796326794895731052901991514519103193L, -3.535050620855721078027883819436759661753e1L);
+  TEST_c_c (cacos, -1.0L, 0x1p50L, 1.570796326794897507409741391764983781004L, -3.535050620855721078027883819436759661753e1L);
+  TEST_c_c (cacos, 1.0L, -0x1p50L, 1.570796326794895731052901991514519103193L, 3.535050620855721078027883819436759661753e1L);
+  TEST_c_c (cacos, -1.0L, -0x1p50L, 1.570796326794897507409741391764983781004L, 3.535050620855721078027883819436759661753e1L);
+#ifndef TEST_FLOAT
+  TEST_c_c (cacos, 0x1p500L, 1.0L, 3.054936363499604682051979393213617699789e-151L, -3.472667374605326000180332928505464606058e2L);
+  TEST_c_c (cacos, 0x1p500L, -1.0L, 3.054936363499604682051979393213617699789e-151L, 3.472667374605326000180332928505464606058e2L);
+  TEST_c_c (cacos, -0x1p500L, 1.0L, 3.141592653589793238462643383279502884197L, -3.472667374605326000180332928505464606058e2L);
+  TEST_c_c (cacos, -0x1p500L, -1.0L, 3.141592653589793238462643383279502884197L, 3.472667374605326000180332928505464606058e2L);
+  TEST_c_c (cacos, 1.0L, 0x1p500L, 1.570796326794896619231321691639751442099L, -3.472667374605326000180332928505464606058e2L);
+  TEST_c_c (cacos, -1.0L, 0x1p500L, 1.570796326794896619231321691639751442099L, -3.472667374605326000180332928505464606058e2L);
+  TEST_c_c (cacos, 1.0L, -0x1p500L, 1.570796326794896619231321691639751442099L, 3.472667374605326000180332928505464606058e2L);
+  TEST_c_c (cacos, -1.0L, -0x1p500L, 1.570796326794896619231321691639751442099L, 3.472667374605326000180332928505464606058e2L);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (cacos, 0x1p5000L, 1.0L, 7.079811261048172892385615158694057552948e-1506L, -3.466429049980286492395577839412341016946e3L);
+  TEST_c_c (cacos, 0x1p5000L, -1.0L, 7.079811261048172892385615158694057552948e-1506L, 3.466429049980286492395577839412341016946e3L);
+  TEST_c_c (cacos, -0x1p5000L, 1.0L, 3.141592653589793238462643383279502884197L, -3.466429049980286492395577839412341016946e3L);
+  TEST_c_c (cacos, -0x1p5000L, -1.0L, 3.141592653589793238462643383279502884197L, 3.466429049980286492395577839412341016946e3L);
+  TEST_c_c (cacos, 1.0L, 0x1p5000L, 1.570796326794896619231321691639751442099L, -3.466429049980286492395577839412341016946e3L);
+  TEST_c_c (cacos, -1.0L, 0x1p5000L, 1.570796326794896619231321691639751442099L, -3.466429049980286492395577839412341016946e3L);
+  TEST_c_c (cacos, 1.0L, -0x1p5000L, 1.570796326794896619231321691639751442099L, 3.466429049980286492395577839412341016946e3L);
+  TEST_c_c (cacos, -1.0L, -0x1p5000L, 1.570796326794896619231321691639751442099L, 3.466429049980286492395577839412341016946e3L);
+#endif
+
+  TEST_c_c (cacos, 0x1.fp127L, 0x1.fp127L, 7.853981633974483096156608458198757210493e-1L, -8.973081118419833726837456344608533993585e1L);
+#ifndef TEST_FLOAT
+  TEST_c_c (cacos, 0x1.fp1023L, 0x1.fp1023L, 7.853981633974483096156608458198757210493e-1L, -7.107906849659093345062145442726115449315e2L);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (cacos, 0x1.fp16383L, 0x1.fp16383L, 7.853981633974483096156608458198757210493e-1L, -1.135753137836666928715489992987020363057e4L);
+#endif
+
   TEST_c_c (cacos, 0.75L, 1.25L, 1.11752014915610270578240049553777969L, -1.13239363160530819522266333696834467L);
   TEST_c_c (cacos, -2, -3, 2.1414491111159960199416055713254211L, 1.9833870299165354323470769028940395L);
 
diff --git a/math/s_cacos.c b/math/s_cacos.c
index 6604b5aec6..acd9b2462a 100644
--- a/math/s_cacos.c
+++ b/math/s_cacos.c
@@ -25,11 +25,27 @@ __cacos (__complex__ double x)
 {
   __complex__ double y;
   __complex__ double res;
-
-  y = __casin (x);
-
-  __real__ res = (double) M_PI_2 - __real__ y;
-  __imag__ res = -__imag__ y;
+  int rcls = fpclassify (__real__ x);
+  int icls = fpclassify (__imag__ x);
+
+  if (rcls <= FP_INFINITE || icls <= FP_INFINITE
+      || (rcls == FP_ZERO && icls == FP_ZERO))
+    {
+      y = __casin (x);
+
+      __real__ res = (double) M_PI_2 - __real__ y;
+      __imag__ res = -__imag__ y;
+    }
+  else
+    {
+      __real__ y = -__imag__ x;
+      __imag__ y = __real__ x;
+
+      y = __kernel_casinh (y, 1);
+
+      __real__ res = __imag__ y;
+      __imag__ res = __real__ y;
+    }
 
   return res;
 }
diff --git a/math/s_cacosf.c b/math/s_cacosf.c
index 04c13e4fa5..df2bf218a3 100644
--- a/math/s_cacosf.c
+++ b/math/s_cacosf.c
@@ -25,11 +25,27 @@ __cacosf (__complex__ float x)
 {
   __complex__ float y;
   __complex__ float res;
-
-  y = __casinf (x);
-
-  __real__ res = (float) M_PI_2 - __real__ y;
-  __imag__ res = -__imag__ y;
+  int rcls = fpclassify (__real__ x);
+  int icls = fpclassify (__imag__ x);
+
+  if (rcls <= FP_INFINITE || icls <= FP_INFINITE
+      || (rcls == FP_ZERO && icls == FP_ZERO))
+    {
+      y = __casinf (x);
+
+      __real__ res = (float) M_PI_2 - __real__ y;
+      __imag__ res = -__imag__ y;
+    }
+  else
+    {
+      __real__ y = -__imag__ x;
+      __imag__ y = __real__ x;
+
+      y = __kernel_casinhf (y, 1);
+
+      __real__ res = __imag__ y;
+      __imag__ res = __real__ y;
+    }
 
   return res;
 }
diff --git a/math/s_cacosl.c b/math/s_cacosl.c
index 304076ddfe..8eab1f0004 100644
--- a/math/s_cacosl.c
+++ b/math/s_cacosl.c
@@ -25,11 +25,27 @@ __cacosl (__complex__ long double x)
 {
   __complex__ long double y;
   __complex__ long double res;
-
-  y = __casinl (x);
-
-  __real__ res = M_PI_2l - __real__ y;
-  __imag__ res = -__imag__ y;
+  int rcls = fpclassify (__real__ x);
+  int icls = fpclassify (__imag__ x);
+
+  if (rcls <= FP_INFINITE || icls <= FP_INFINITE
+      || (rcls == FP_ZERO && icls == FP_ZERO))
+    {
+      y = __casinl (x);
+
+      __real__ res = M_PI_2l - __real__ y;
+      __imag__ res = -__imag__ y;
+    }
+  else
+    {
+      __real__ y = -__imag__ x;
+      __imag__ y = __real__ x;
+
+      y = __kernel_casinhl (y, 1);
+
+      __real__ res = __imag__ y;
+      __imag__ res = __real__ y;
+    }
 
   return res;
 }
diff --git a/math/s_casinh.c b/math/s_casinh.c
index b493982d88..657e269ac1 100644
--- a/math/s_casinh.c
+++ b/math/s_casinh.c
@@ -20,7 +20,6 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-#include <float.h>
 
 __complex__ double
 __casinh (__complex__ double x)
@@ -62,40 +61,7 @@ __casinh (__complex__ double x)
     }
   else
     {
-      double rx, ix;
-      __complex__ double y;
-
-      /* Avoid cancellation by reducing to the first quadrant.  */
-      rx = fabs (__real__ x);
-      ix = fabs (__imag__ x);
-
-      if (rx >= 1.0 / DBL_EPSILON || ix >= 1.0 / DBL_EPSILON)
-	{
-	  /* For large x in the first quadrant, x + csqrt (1 + x * x)
-	     is sufficiently close to 2 * x to make no significant
-	     difference to the result; avoid possible overflow from
-	     the squaring and addition.  */
-	  __real__ y = rx;
-	  __imag__ y = ix;
-	  res = __clog (y);
-	  __real__ res += M_LN2;
-	}
-      else
-	{
-	  __real__ y = (rx - ix) * (rx + ix) + 1.0;
-	  __imag__ y = 2.0 * rx * ix;
-
-	  y = __csqrt (y);
-
-	  __real__ y += rx;
-	  __imag__ y += ix;
-
-	  res = __clog (y);
-	}
-
-      /* Give results the correct sign for the original argument.  */
-      __real__ res = __copysign (__real__ res, __real__ x);
-      __imag__ res = __copysign (__imag__ res, __imag__ x);
+      res = __kernel_casinh (x, 0);
     }
 
   return res;
diff --git a/math/s_casinhf.c b/math/s_casinhf.c
index f865e14490..8663c2e7cc 100644
--- a/math/s_casinhf.c
+++ b/math/s_casinhf.c
@@ -20,7 +20,6 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-#include <float.h>
 
 __complex__ float
 __casinhf (__complex__ float x)
@@ -62,40 +61,7 @@ __casinhf (__complex__ float x)
     }
   else
     {
-      float rx, ix;
-      __complex__ float y;
-
-      /* Avoid cancellation by reducing to the first quadrant.  */
-      rx = fabsf (__real__ x);
-      ix = fabsf (__imag__ x);
-
-      if (rx >= 1.0f / FLT_EPSILON || ix >= 1.0f / FLT_EPSILON)
-	{
-	  /* For large x in the first quadrant, x + csqrt (1 + x * x)
-	     is sufficiently close to 2 * x to make no significant
-	     difference to the result; avoid possible overflow from
-	     the squaring and addition.  */
-	  __real__ y = rx;
-	  __imag__ y = ix;
-	  res = __clogf (y);
-	  __real__ res += (float) M_LN2;
-	}
-      else
-	{
-	  __real__ y = (rx - ix) * (rx + ix) + 1.0;
-	  __imag__ y = 2.0 * rx * ix;
-
-	  y = __csqrtf (y);
-
-	  __real__ y += rx;
-	  __imag__ y += ix;
-
-	  res = __clogf (y);
-	}
-
-      /* Give results the correct sign for the original argument.  */
-      __real__ res = __copysignf (__real__ res, __real__ x);
-      __imag__ res = __copysignf (__imag__ res, __imag__ x);
+      res = __kernel_casinhf (x, 0);
     }
 
   return res;
diff --git a/math/s_casinhl.c b/math/s_casinhl.c
index d7c74593e4..2afc52714e 100644
--- a/math/s_casinhl.c
+++ b/math/s_casinhl.c
@@ -20,14 +20,6 @@
 #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
 __casinhl (__complex__ long double x)
@@ -69,40 +61,7 @@ __casinhl (__complex__ long double x)
     }
   else
     {
-      long double rx, ix;
-      __complex__ long double y;
-
-      /* Avoid cancellation by reducing to the first quadrant.  */
-      rx = fabsl (__real__ x);
-      ix = fabsl (__imag__ x);
-
-      if (rx >= 1.0L / LDBL_EPSILON || ix >= 1.0L / LDBL_EPSILON)
-	{
-	  /* For large x in the first quadrant, x + csqrt (1 + x * x)
-	     is sufficiently close to 2 * x to make no significant
-	     difference to the result; avoid possible overflow from
-	     the squaring and addition.  */
-	  __real__ y = rx;
-	  __imag__ y = ix;
-	  res = __clogl (y);
-	  __real__ res += M_LN2l;
-	}
-      else
-	{
-	  __real__ y = (rx - ix) * (rx + ix) + 1.0;
-	  __imag__ y = 2.0 * rx * ix;
-
-	  y = __csqrtl (y);
-
-	  __real__ y += rx;
-	  __imag__ y += ix;
-
-	  res = __clogl (y);
-	}
-
-      /* Give results the correct sign for the original argument.  */
-      __real__ res = __copysignl (__real__ res, __real__ x);
-      __imag__ res = __copysignl (__imag__ res, __imag__ x);
+      res = __kernel_casinhl (x, 0);
     }
 
   return res;
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 3fc30de462..1525b16f37 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -303,6 +303,12 @@ float: 1
 ifloat: 1
 ildouble: 2
 ldouble: 2
+Test "Imaginary part of: cacos (0x1.fp1023 + 0x1.fp1023 i) == 7.853981633974483096156608458198757210493e-1 - 7.107906849659093345062145442726115449315e2 i":
+double: 1
+idouble: 1
+Test "Imaginary part of: cacos (0x1.fp127 + 0x1.fp127 i) == 7.853981633974483096156608458198757210493e-1 - 8.973081118419833726837456344608533993585e1 i":
+double: 1
+idouble: 1
 Test "Imaginary part of: cacos (1.5 + +0 i) == +0 - 0.9624236501192068949955178268487368462704 i":
 double: 1
 float: 1
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 95b6aec81e..63c6aed2a5 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -244,6 +244,12 @@ ifloat: 1
 Test "Imaginary part of: cacos (-0 - 1.5 i) == pi/2 + 1.194763217287109304111930828519090523536 i":
 double: 1
 idouble: 1
+Test "Real part of: cacos (-1.0 + 0x1p50 i) == 1.570796326794897507409741391764983781004 - 3.535050620855721078027883819436759661753e1 i":
+float: 1
+ifloat: 1
+Test "Real part of: cacos (-1.0 - 0x1p50 i) == 1.570796326794897507409741391764983781004 + 3.535050620855721078027883819436759661753e1 i":
+float: 1
+ifloat: 1
 Test "Imaginary part of: cacos (-1.5 + +0 i) == pi - 0.9624236501192068949955178268487368462704 i":
 double: 1
 float: 1
@@ -254,6 +260,9 @@ ldouble: 1
 Test "Imaginary part of: cacos (-1.5 - 0 i) == pi + 0.9624236501192068949955178268487368462704 i":
 ildouble: 1
 ldouble: 1
+Test "Real part of: cacos (-2 - 3 i) == 2.1414491111159960199416055713254211 + 1.9833870299165354323470769028940395 i":
+float: 1
+ifloat: 1
 Test "Real part of: cacos (0.5 + +0 i) == 1.047197551196597746154214461093167628066 - 0 i":
 double: 1
 idouble: 1
@@ -272,6 +281,12 @@ float: 1
 ifloat: 1
 ildouble: 2
 ldouble: 2
+Test "Imaginary part of: cacos (0x1.fp1023 + 0x1.fp1023 i) == 7.853981633974483096156608458198757210493e-1 - 7.107906849659093345062145442726115449315e2 i":
+double: 1
+idouble: 1
+Test "Imaginary part of: cacos (0x1.fp127 + 0x1.fp127 i) == 7.853981633974483096156608458198757210493e-1 - 8.973081118419833726837456344608533993585e1 i":
+double: 1
+idouble: 1
 Test "Imaginary part of: cacos (1.5 + +0 i) == +0 - 0.9624236501192068949955178268487368462704 i":
 double: 1
 float: 1