Fix cacos real-part inaccuracy for result real part near 0 (bug 15023).
authorJoseph Myers <joseph@codesourcery.com>
Thu, 17 Jan 2013 20:25:51 +0000 (20:25 +0000)
committerJoseph Myers <joseph@codesourcery.com>
Thu, 17 Jan 2013 20:25:51 +0000 (20:25 +0000)
16 files changed:
ChangeLog
NEWS
include/complex.h
math/Makefile
math/k_casinh.c [new file with mode: 0644]
math/k_casinhf.c [new file with mode: 0644]
math/k_casinhl.c [new file with mode: 0644]
math/libm-test.inc
math/s_cacos.c
math/s_cacosf.c
math/s_cacosl.c
math/s_casinh.c
math/s_casinhf.c
math/s_casinhl.c
sysdeps/i386/fpu/libm-test-ulps
sysdeps/x86_64/fpu/libm-test-ulps

index e7cfb64..890e3d4 100644 (file)
--- 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 55f65fd..99fb542 100644 (file)
--- 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.
 
 \f
 Version 2.17
index acf8cf1..e173f1f 100644 (file)
@@ -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 b9519cf..da18b56 100644 (file)
@@ -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 (file)
index 0000000..7f98f24
--- /dev/null
@@ -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 (file)
index 0000000..9401636
--- /dev/null
@@ -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 (file)
index 0000000..6412979
--- /dev/null
@@ -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 56e3217..1c2970f 100644 (file)
@@ -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);
 
index 6604b5a..acd9b24 100644 (file)
@@ -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;
 }
index 04c13e4..df2bf21 100644 (file)
@@ -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;
 }
index 304076d..8eab1f0 100644 (file)
@@ -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;
 }
index b493982..657e269 100644 (file)
@@ -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;
index f865e14..8663c2e 100644 (file)
@@ -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;
index d7c7459..2afc527 100644 (file)
 #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;
index 3fc30de..1525b16 100644 (file)
@@ -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
index 95b6aec..63c6aed 100644 (file)
@@ -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