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.

Index: glibc-2.11.3/include/complex.h
===================================================================
--- glibc-2.11.3.orig/include/complex.h
+++ glibc-2.11.3/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
Index: glibc-2.11.3/math/Makefile
===================================================================
--- glibc-2.11.3.orig/math/Makefile
+++ glibc-2.11.3/math/Makefile
@@ -60,7 +60,7 @@ libm-calls = e_acos e_acosh e_asin e_ata
 	     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	\
-	     $(calls:s_%=m_%)
+	     $(calls:s_%=m_%) k_casinh
 
 include ../Makeconfig
 
Index: glibc-2.11.3/math/k_casinh.c
===================================================================
--- /dev/null
+++ glibc-2.11.3/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;
+}
Index: glibc-2.11.3/math/k_casinhf.c
===================================================================
--- /dev/null
+++ glibc-2.11.3/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;
+}
Index: glibc-2.11.3/math/k_casinhl.c
===================================================================
--- /dev/null
+++ glibc-2.11.3/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;
+}
Index: glibc-2.11.3/math/libm-test.inc
===================================================================
--- glibc-2.11.3.orig/math/libm-test.inc
+++ glibc-2.11.3/math/libm-test.inc
@@ -1092,6 +1092,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);
 
Index: glibc-2.11.3/math/s_cacos.c
===================================================================
--- glibc-2.11.3.orig/math/s_cacos.c
+++ glibc-2.11.3/math/s_cacos.c
@@ -26,11 +26,27 @@ __cacos (__complex__ double x)
 {
   __complex__ double y;
   __complex__ double res;
+  int rcls = fpclassify (__real__ x);
+  int icls = fpclassify (__imag__ x);
 
-  y = __casin (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;
+      __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;
 }
Index: glibc-2.11.3/math/s_cacosf.c
===================================================================
--- glibc-2.11.3.orig/math/s_cacosf.c
+++ glibc-2.11.3/math/s_cacosf.c
@@ -26,11 +26,27 @@ __cacosf (__complex__ float x)
 {
   __complex__ float y;
   __complex__ float res;
+  int rcls = fpclassify (__real__ x);
+  int icls = fpclassify (__imag__ x);
 
-  y = __casinf (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;
+      __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;
 }
Index: glibc-2.11.3/math/s_cacosl.c
===================================================================
--- glibc-2.11.3.orig/math/s_cacosl.c
+++ glibc-2.11.3/math/s_cacosl.c
@@ -26,11 +26,27 @@ __cacosl (__complex__ long double x)
 {
   __complex__ long double y;
   __complex__ long double res;
+  int rcls = fpclassify (__real__ x);
+  int icls = fpclassify (__imag__ x);
 
-  y = __casinl (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;
+      __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;
 }
Index: glibc-2.11.3/math/s_casinh.c
===================================================================
--- glibc-2.11.3.orig/math/s_casinh.c
+++ glibc-2.11.3/math/s_casinh.c
@@ -20,7 +20,6 @@
 
 #include <complex.h>
 #include <math.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;
Index: glibc-2.11.3/math/s_casinhf.c
===================================================================
--- glibc-2.11.3.orig/math/s_casinhf.c
+++ glibc-2.11.3/math/s_casinhf.c
@@ -20,7 +20,6 @@
 
 #include <complex.h>
 #include <math.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;
Index: glibc-2.11.3/math/s_casinhl.c
===================================================================
--- glibc-2.11.3.orig/math/s_casinhl.c
+++ glibc-2.11.3/math/s_casinhl.c
@@ -20,14 +20,6 @@
 
 #include <complex.h>
 #include <math.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;
Index: glibc-2.11.3/sysdeps/i386/fpu/libm-test-ulps
===================================================================
--- glibc-2.11.3.orig/sysdeps/i386/fpu/libm-test-ulps
+++ glibc-2.11.3/sysdeps/i386/fpu/libm-test-ulps
@@ -116,6 +116,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
Index: glibc-2.11.3/sysdeps/x86_64/fpu/libm-test-ulps
===================================================================
--- glibc-2.11.3.orig/sysdeps/x86_64/fpu/libm-test-ulps
+++ glibc-2.11.3/sysdeps/x86_64/fpu/libm-test-ulps
@@ -101,6 +101,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
@@ -111,6 +117,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
@@ -129,6 +138,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
