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

11 files changed, 366 insertions, 128 deletions

diff --git a/math/Makefile b/math/Makefile
index b9519cf..da18b56 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 0000000..7f98f24
--- /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 0000000..9401636
--- /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 0000000..6412979
--- /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 56e3217..1c2970f 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 6604b5a..acd9b24 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 04c13e4..df2bf21 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 304076d..8eab1f0 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 b493982..657e269 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 f865e14..8663c2e 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 d7c7459..2afc527 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;