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

7 files changed, 109 insertions, 9 deletions

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;