This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix catan, catanh inaccuracy through use of log (bug 15394)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Wed, 24 Apr 2013 16:42:01 +0000
- Subject: Fix catan, catanh inaccuracy through use of log (bug 15394)
Bug 15394 is inaccuracy of catan and catanh for certain arguments
where they involve computing log of a number close to 1 (or in the
case of catanh, subtracting logs of similar arguments, which could
lead to worse results in some cases than computing the log of their
ratio, through large cancellation even when the ratio isn't that close
to 1).
This patch fixes this inaccuracy by using log1p except in the cases
where log of a number close to 0 is involved, when use of log1p of a
number close to -1 would yield worse results than log and so log of a
ratio is used instead.
Tested x86_64 and x86 and ulps updated accordingly.
Note that there are plenty more issues with catan/catanh that I'll
address separately.
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.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 0049fcd..447b603 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 46c18bf..783941a 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 5a43247..0dc85ff 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 5371f44..0ee8c64 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 8385af4..ca9a301 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 6844f03..4897c0c 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 57d2e59..e04dba7 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 fecaa94..a6d9b0f 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 9c9c473..2747fc5 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
--
Joseph S. Myers
joseph@codesourcery.com