Fix catan, catanh inaccuracy through use of log (bug 15394).

11 files changed, 130 insertions, 10 deletions

diff --git a/ChangeLog b/ChangeLog
index 2d44b8ed60..c3404acd82 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,19 @@
+2013-04-24  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #15394]
+	* math/s_catan.c (__catan): Calculate imaginary part of result
+	with log1p not log unless computing log of number close to 0.
+	* math/s_catanf.c (__catanf): Likewise.
+	* math/s_catanl.c (__catanl): Likewise.
+	* math/s_catanh.c (__catanh): Calculate real part of result with
+	log1p not log unless computing log of number close to 0.
+	* math/s_catanhf.c (__catanhf): Likewise.
+	* math/s_catanhl.c (__catanhl): Likewise.
+	* math/libm-test.inc (catan_test): Add more tests.
+	(catanh_test): Likewise.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+
 2013-04-24  Siddhesh Poyarekar  <siddhesh@redhat.com>
 
 	* benchtests/Makefile: Mention files in which fast and slow
diff --git a/NEWS b/NEWS
index a1aeb42550..e5c4f908ef 100644
--- a/NEWS
+++ b/NEWS
@@ -14,7 +14,7 @@ Version 2.18
   14888, 14920, 14964, 14981, 14982, 14985, 14994, 14996, 15003, 15006,
   15020, 15023, 15036, 15054, 15055, 15062, 15078, 15160, 15214, 15232,
   15234, 15283, 15285, 15287, 15304, 15305, 15307, 15309, 15327, 15330,
-  15335, 15336, 15337, 15342, 15346, 15361.
+  15335, 15336, 15337, 15342, 15346, 15361, 15394.
 
 * CVE-2013-0242 Buffer overrun in regexp matcher has been fixed (Bugzilla
   #15078).
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 0049fcda38..447b6037ae 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -4299,6 +4299,35 @@ catan_test (void)
 
   TEST_c_c (catan, qnan_value, qnan_value, qnan_value, qnan_value);
 
+  TEST_c_c (catan, 0x1p50L, 0.0L, 1.570796326794895731052901991514519103193L, 0.0L);
+  TEST_c_c (catan, 0x1p50L, -0.0L, 1.570796326794895731052901991514519103193L, -0.0L);
+  TEST_c_c (catan, -0x1p50L, 0.0L, -1.570796326794895731052901991514519103193L, 0.0L);
+  TEST_c_c (catan, -0x1p50L, -0.0L, -1.570796326794895731052901991514519103193L, -0.0L);
+  TEST_c_c (catan, 0.0L, 0x1p50L, 1.570796326794896619231321691639751442099L, 8.881784197001252323389053344728897997441e-16L);
+  TEST_c_c (catan, -0.0L, 0x1p50L, -1.570796326794896619231321691639751442099L, 8.881784197001252323389053344728897997441e-16L);
+  TEST_c_c (catan, 0.0L, -0x1p50L, 1.570796326794896619231321691639751442099L, -8.881784197001252323389053344728897997441e-16L);
+  TEST_c_c (catan, -0.0L, -0x1p50L, -1.570796326794896619231321691639751442099L, -8.881784197001252323389053344728897997441e-16L);
+#ifndef TEST_FLOAT
+  TEST_c_c (catan, 0x1p500L, 0.0L, 1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, 0x1p500L, -0.0L, 1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, -0x1p500L, 0.0L, -1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, -0x1p500L, -0.0L, -1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, 0.0L, 0x1p500L, 1.570796326794896619231321691639751442099L, 3.054936363499604682051979393213617699789e-151L);
+  TEST_c_c (catan, -0.0L, 0x1p500L, -1.570796326794896619231321691639751442099L, 3.054936363499604682051979393213617699789e-151L);
+  TEST_c_c (catan, 0.0L, -0x1p500L, 1.570796326794896619231321691639751442099L, -3.054936363499604682051979393213617699789e-151L);
+  TEST_c_c (catan, -0.0L, -0x1p500L, -1.570796326794896619231321691639751442099L, -3.054936363499604682051979393213617699789e-151L);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (catan, 0x1p5000L, 0.0L, 1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, 0x1p5000L, -0.0L, 1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, -0x1p5000L, 0.0L, -1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, -0x1p5000L, -0.0L, -1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, 0.0L, 0x1p5000L, 1.570796326794896619231321691639751442099L, 7.079811261048172892385615158694057552948e-1506L);
+  TEST_c_c (catan, -0.0L, 0x1p5000L, -1.570796326794896619231321691639751442099L, 7.079811261048172892385615158694057552948e-1506L);
+  TEST_c_c (catan, 0.0L, -0x1p5000L, 1.570796326794896619231321691639751442099L, -7.079811261048172892385615158694057552948e-1506L);
+  TEST_c_c (catan, -0.0L, -0x1p5000L, -1.570796326794896619231321691639751442099L, -7.079811261048172892385615158694057552948e-1506L);
+#endif
+
   TEST_c_c (catan, 0.75L, 1.25L, 1.10714871779409050301706546017853704L, 0.549306144334054845697622618461262852L);
   TEST_c_c (catan, -2, -3, -1.4099210495965755225306193844604208L, -0.22907268296853876629588180294200276L);
 
@@ -4365,6 +4394,35 @@ catanh_test (void)
 
   TEST_c_c (catanh, qnan_value, qnan_value, qnan_value, qnan_value);
 
+  TEST_c_c (catanh, 0x1p50L, 0.0L, 8.881784197001252323389053344728897997441e-16L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1p50L, -0.0L, 8.881784197001252323389053344728897997441e-16L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p50L, 0.0L, -8.881784197001252323389053344728897997441e-16L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p50L, -0.0L, -8.881784197001252323389053344728897997441e-16L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, 0x1p50L, 0.0L, 1.570796326794895731052901991514519103193L);
+  TEST_c_c (catanh, -0.0L, 0x1p50L, -0.0L, 1.570796326794895731052901991514519103193L);
+  TEST_c_c (catanh, 0.0L, -0x1p50L, 0.0L, -1.570796326794895731052901991514519103193L);
+  TEST_c_c (catanh, -0.0L, -0x1p50L, -0.0L, -1.570796326794895731052901991514519103193L);
+#ifndef TEST_FLOAT
+  TEST_c_c (catanh, 0x1p500L, 0.0L, 3.054936363499604682051979393213617699789e-151L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1p500L, -0.0L, 3.054936363499604682051979393213617699789e-151L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p500L, 0.0L, -3.054936363499604682051979393213617699789e-151L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p500L, -0.0L, -3.054936363499604682051979393213617699789e-151L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, 0x1p500L, 0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, 0x1p500L, -0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, -0x1p500L, 0.0L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, -0x1p500L, -0.0L, -1.570796326794896619231321691639751442099L);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (catanh, 0x1p5000L, 0.0L, 7.079811261048172892385615158694057552948e-1506L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1p5000L, -0.0L, 7.079811261048172892385615158694057552948e-1506L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p5000L, 0.0L, -7.079811261048172892385615158694057552948e-1506L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p5000L, -0.0L, -7.079811261048172892385615158694057552948e-1506L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, 0x1p5000L, 0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, 0x1p5000L, -0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, -0x1p5000L, 0.0L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, -0x1p5000L, -0.0L, -1.570796326794896619231321691639751442099L);
+#endif
+
   TEST_c_c (catanh, 0.75L, 1.25L, 0.261492138795671927078652057366532140L, 0.996825126463918666098902241310446708L);
   TEST_c_c (catanh, -2, -3, -0.14694666622552975204743278515471595L, -1.3389725222944935611241935759091443L);
 
diff --git a/math/s_catan.c b/math/s_catan.c
index 46c18bf2af..783941a72e 100644
--- a/math/s_catan.c
+++ b/math/s_catan.c
@@ -61,7 +61,7 @@ __catan (__complex__ double x)
     }
   else
     {
-      double r2, num, den;
+      double r2, num, den, f;
 
       r2 = __real__ x * __real__ x;
 
@@ -75,7 +75,14 @@ __catan (__complex__ double x)
       den = __imag__ x - 1.0;
       den = r2 + den * den;
 
-      __imag__ res = 0.25 * __ieee754_log (num / den);
+      f = num / den;
+      if (f < 0.5)
+	__imag__ res = 0.25 * __ieee754_log (f);
+      else
+	{
+	  num = 4.0 * __imag__ x;
+	  __imag__ res = 0.25 * __log1p (num / den);
+	}
     }
 
   return res;
diff --git a/math/s_catanf.c b/math/s_catanf.c
index 5a432471cd..0dc85ffc25 100644
--- a/math/s_catanf.c
+++ b/math/s_catanf.c
@@ -61,7 +61,7 @@ __catanf (__complex__ float x)
     }
   else
     {
-      float r2, num, den;
+      float r2, num, den, f;
 
       r2 = __real__ x * __real__ x;
 
@@ -75,7 +75,14 @@ __catanf (__complex__ float x)
       den = __imag__ x - 1.0;
       den = r2 + den * den;
 
-      __imag__ res = 0.25 * __ieee754_logf (num / den);
+      f = num / den;
+      if (f < 0.5)
+	__imag__ res = 0.25 * __ieee754_logf (f);
+      else
+	{
+	  num = 4.0 * __imag__ x;
+	  __imag__ res = 0.25 * __log1pf (num / den);
+	}
     }
 
   return res;
diff --git a/math/s_catanh.c b/math/s_catanh.c
index 5371f44b67..0ee8c64b6f 100644
--- a/math/s_catanh.c
+++ b/math/s_catanh.c
@@ -64,7 +64,14 @@ __catanh (__complex__ double x)
       double den = 1.0 - __real__ x;
       den = i2 + den * den;
 
-      __real__ res = 0.25 * (__ieee754_log (num) - __ieee754_log (den));
+      double f = num / den;
+      if (f < 0.5)
+	__real__ res = 0.25 * __ieee754_log (f);
+      else
+	{
+	  num = 4.0 * __real__ x;
+	  __real__ res = 0.25 * __log1p (num / den);
+	}
 
       den = 1 - __real__ x * __real__ x - i2;
 
diff --git a/math/s_catanhf.c b/math/s_catanhf.c
index 8385af41b0..ca9a30101e 100644
--- a/math/s_catanhf.c
+++ b/math/s_catanhf.c
@@ -64,7 +64,14 @@ __catanhf (__complex__ float x)
       float den = 1.0 - __real__ x;
       den = i2 + den * den;
 
-      __real__ res = 0.25 * (__ieee754_logf (num) - __ieee754_logf (den));
+      float f = num / den;
+      if (f < 0.5)
+	__real__ res = 0.25 * __ieee754_logf (f);
+      else
+	{
+	  num = 4.0 * __real__ x;
+	  __real__ res = 0.25 * __log1pf (num / den);
+	}
 
       den = 1 - __real__ x * __real__ x - i2;
 
diff --git a/math/s_catanhl.c b/math/s_catanhl.c
index 6844f03dea..4897c0c871 100644
--- a/math/s_catanhl.c
+++ b/math/s_catanhl.c
@@ -64,7 +64,14 @@ __catanhl (__complex__ long double x)
       long double den = 1.0 - __real__ x;
       den = i2 + den * den;
 
-      __real__ res = 0.25 * (__ieee754_logl (num) - __ieee754_logl (den));
+      long double f = num / den;
+      if (f < 0.5)
+	__real__ res = 0.25 * __ieee754_logl (f);
+      else
+	{
+	  num = 4.0 * __real__ x;
+	  __real__ res = 0.25 * __log1pl (num / den);
+	}
 
       den = 1 - __real__ x * __real__ x - i2;
 
diff --git a/math/s_catanl.c b/math/s_catanl.c
index 57d2e59550..e04dba7387 100644
--- a/math/s_catanl.c
+++ b/math/s_catanl.c
@@ -61,7 +61,7 @@ __catanl (__complex__ long double x)
     }
   else
     {
-      long double r2, num, den;
+      long double r2, num, den, f;
 
       r2 = __real__ x * __real__ x;
 
@@ -75,7 +75,14 @@ __catanl (__complex__ long double x)
       den = __imag__ x - 1.0;
       den = r2 + den * den;
 
-      __imag__ res = 0.25 * __ieee754_logl (num / den);
+      f = num / den;
+      if (f < 0.5)
+	__imag__ res = 0.25 * __ieee754_logl (f);
+      else
+	{
+	  num = 4.0 * __imag__ x;
+	  __imag__ res = 0.25 * __log1pl (num / den);
+	}
     }
 
   return res;
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index fecaa94d45..a6d9b0f247 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -3417,6 +3417,8 @@ ldouble: 1
 Test "Real part of: catanh (0.75 + 1.25 i) == 0.261492138795671927078652057366532140 + 0.996825126463918666098902241310446708 i":
 double: 1
 idouble: 1
+ildouble: 1
+ldouble: 1
 
 # cbrt
 Test "cbrt (-27.0) == -3.0":
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 9c9c473449..2747fc51cd 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -3985,6 +3985,8 @@ ifloat: 4
 Test "Real part of: catanh (0.75 + 1.25 i) == 0.261492138795671927078652057366532140 + 0.996825126463918666098902241310446708 i":
 double: 1
 idouble: 1
+ildouble: 1
+ldouble: 1
 Test "Imaginary part of: catanh (0.75 + 1.25 i) == 0.261492138795671927078652057366532140 + 0.996825126463918666098902241310446708 i":
 float: 6
 ifloat: 6