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

ChangeLog

@@ -1,3 +1,18 @@
+2013-03-30  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #10357]
+	* math/k_casinh.c (__kernel_casinh): Handle arguments with
+	imaginary part less than 1.0 and real part less than 0.5
+	specially.
+	* math/k_casinhf.c (__kernel_casinhf): Likewise.
+	* math/k_casinhl.c (__kernel_casinhl): Likewise.
+	* math/libm-test.inc (UNDERFLOW_EXCEPTION_OK_DOUBLE): New macro.
+	(cacos_test): Add more tests.
+	(casin_test): Likewise.
+	(casinh_test): Likewise.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+
 2013-03-29  Siddhesh Poyarekar  <siddhesh@redhat.com>
 
 	* sysdeps/powerpc/power4/fpu/mpa-arch.h (INTEGER_OF): Replace

NEWS

@@ -9,10 +9,10 @@ Version 2.18
 
 * The following bugs are resolved with this release:
 
-  11120, 11561, 12723, 13550, 13889, 13951, 14142, 14176, 14200, 14317,
+  10357, 11120, 11561, 12723, 13550, 13889, 13951, 14142, 14176, 14200,
-  14327, 14496, 14812, 14920, 14964, 14981, 14982, 14985, 14994, 14996,
+  14317, 14327, 14496, 14812, 14920, 14964, 14981, 14982, 14985, 14994,
-  15003, 15006, 15020, 15023, 15036, 15054, 15055, 15062, 15078, 15160,
+  14996, 15003, 15006, 15020, 15023, 15036, 15054, 15055, 15062, 15078,
-  15214, 15232, 15234, 15283, 15285, 15287, 15304, 15307.
+  15160, 15214, 15232, 15234, 15283, 15285, 15287, 15304, 15307.
 
 * Add support for calling C++11 thread_local object destructors on thread
   and program exit.  This needs compiler support for offloading C++11

math/k_casinh.c

@@ -134,6 +134,58 @@ __kernel_casinh (__complex__ double x, int adj)
 	    __imag__ res = __ieee754_atan2 (1.0 + s2, rx + s1);
 	}
     }
+  else if (ix < 1.0 && rx < 0.5)
+    {
+      if (ix >= DBL_EPSILON)
+	{
+	  if (rx < DBL_EPSILON * DBL_EPSILON)
+	    {
+	      double onemix2 = (1.0 + ix) * (1.0 - ix);
+	      double s = __ieee754_sqrt (onemix2);
+
+	      __real__ res = __log1p (2.0 * rx / s) / 2.0;
+	      if (adj)
+		__imag__ res = __ieee754_atan2 (s, __imag__ x);
+	      else
+		__imag__ res = __ieee754_atan2 (ix, s);
+	    }
+	  else
+	    {
+	      double onemix2 = (1.0 + ix) * (1.0 - ix);
+	      double rx2 = rx * rx;
+	      double f = rx2 * (2.0 + rx2 + 2.0 * ix * ix);
+	      double d = __ieee754_sqrt (onemix2 * onemix2 + f);
+	      double dp = d + onemix2;
+	      double dm = f / dp;
+	      double r1 = __ieee754_sqrt ((dp + rx2) / 2.0);
+	      double r2 = rx * ix / r1;
+
+	      __real__ res
+		= __log1p (rx2 + dm + 2.0 * (rx * r1 + ix * r2)) / 2.0;
+	      if (adj)
+		__imag__ res = __ieee754_atan2 (rx + r1,
+						__copysign (ix + r2,
+							    __imag__ x));
+	      else
+		__imag__ res = __ieee754_atan2 (ix + r2, rx + r1);
+	    }
+	}
+      else
+	{
+	  double s = __ieee754_hypot (1.0, rx);
+
+	  __real__ res = __log1p (2.0 * rx * (rx + s)) / 2.0;
+	  if (adj)
+	    __imag__ res = __ieee754_atan2 (s, __imag__ x);
+	  else
+	    __imag__ res = __ieee754_atan2 (ix, s);
+	}
+      if (__real__ res < DBL_MIN)
+	{
+	  volatile double force_underflow = __real__ res * __real__ res;
+	  (void) force_underflow;
+	}
+    }
   else
     {
       __real__ y = (rx - ix) * (rx + ix) + 1.0;

math/k_casinhf.c

@@ -136,6 +136,58 @@ __kernel_casinhf (__complex__ float x, int adj)
 	    __imag__ res = __ieee754_atan2f (1.0f + s2, rx + s1);
 	}
     }
+  else if (ix < 1.0f && rx < 0.5f)
+    {
+      if (ix >= FLT_EPSILON)
+	{
+	  if (rx < FLT_EPSILON * FLT_EPSILON)
+	    {
+	      float onemix2 = (1.0f + ix) * (1.0f - ix);
+	      float s = __ieee754_sqrtf (onemix2);
+
+	      __real__ res = __log1pf (2.0f * rx / s) / 2.0f;
+	      if (adj)
+		__imag__ res = __ieee754_atan2f (s, __imag__ x);
+	      else
+		__imag__ res = __ieee754_atan2f (ix, s);
+	    }
+	  else
+	    {
+	      float onemix2 = (1.0f + ix) * (1.0f - ix);
+	      float rx2 = rx * rx;
+	      float f = rx2 * (2.0f + rx2 + 2.0f * ix * ix);
+	      float d = __ieee754_sqrtf (onemix2 * onemix2 + f);
+	      float dp = d + onemix2;
+	      float dm = f / dp;
+	      float r1 = __ieee754_sqrtf ((dp + rx2) / 2.0f);
+	      float r2 = rx * ix / r1;
+
+	      __real__ res
+		= __log1pf (rx2 + dm + 2.0f * (rx * r1 + ix * r2)) / 2.0f;
+	      if (adj)
+		__imag__ res = __ieee754_atan2f (rx + r1,
+						 __copysignf (ix + r2,
+							      __imag__ x));
+	      else
+		__imag__ res = __ieee754_atan2f (ix + r2, rx + r1);
+	    }
+	}
+      else
+	{
+	  float s = __ieee754_hypotf (1.0f, rx);
+
+	  __real__ res = __log1pf (2.0f * rx * (rx + s)) / 2.0f;
+	  if (adj)
+	    __imag__ res = __ieee754_atan2f (s, __imag__ x);
+	  else
+	    __imag__ res = __ieee754_atan2f (ix, s);
+	}
+      if (__real__ res < FLT_MIN)
+	{
+	  volatile float force_underflow = __real__ res * __real__ res;
+	  (void) force_underflow;
+	}
+    }
   else
     {
       __real__ y = (rx - ix) * (rx + ix) + 1.0f;

math/k_casinhl.c

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