This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: Fix ctan, ctanh of subnormals in round-upwards mode (bug 14328)
- From: Adhemerval Zanella <azanella at linux dot vnet dot ibm dot com>
- To: libc-alpha at sourceware dot org
- Date: Sun, 08 Jul 2012 09:57:04 -0300
- Subject: Re: Fix ctan, ctanh of subnormals in round-upwards mode (bug 14328)
- References: <Pine.LNX.4.64.1207032130450.9313@digraph.polyomino.org.uk>
On 07/03/2012 06:32 PM, Joseph S. Myers wrote:
> Bug 14328 is inaccuracy from ctan and ctanh on certain cases of
> subnormal inputs, in round-upwards mode only.
>
This patch is the same correction, but for IBM long double plus PPC
ulps update. The correction relies on LDBL_EPSILON value, however
for IBM long double, GCC builtin sets LDBL_EPSILON as LDBL_DENORM_MIN
and thus the comparison:
if (fabsl (sinhrx) > fabsl (cosix) * LDBL_EPSILON)
Issues a FE_UNDERFLOW exception in the multiplication for some values.
To fix this I used the constant 2^-106 as IBM long double epsilon
(since IBM long double format target this precision). The ULPs for
different rounding are quite high mainly because the IBM long double
arithmetic operators only guarantees good precision in FE_TONEAREST
rounding mode.
Any thoughts, tips, advices?
--
2012-07-08 Adhemerval Zanella <azanella@linux.vnet.ibm.com>
* sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c: Do not call sinh and cosh
for subnormals or multiply small sinh result by itself.
* sysdeps/ieee754/ldbl-128ibm/s_ctanl.c: Likewise.
* sysdeps/powerpc/fpu/libm-test-ulps: Update.
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c b/sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c
index 2ab80a2..43ed6d3 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c
@@ -25,6 +25,8 @@
#include <math_private.h>
+/* IBM long double GCC builtin sets LDBL_EPSILON == LDBL_DENORM_MIN */
+static const long double ldbl_eps = 1.232595164407831e-32; /* 2^-106 */
__complex__ long double
__ctanhl (__complex__ long double x)
@@ -35,8 +37,8 @@ __ctanhl (__complex__ long double x)
{
if (__isinfl (__real__ x))
{
- __real__ res = __copysignl (1.0, __real__ x);
- __imag__ res = __copysignl (0.0, __imag__ x);
+ __real__ res = __copysignl (1.0L, __real__ x);
+ __imag__ res = __copysignl (0.0L, __imag__ x);
}
else if (__imag__ x == 0.0)
{
@@ -57,7 +59,7 @@ __ctanhl (__complex__ long double x)
{
long double sinix, cosix;
long double den;
- const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2);
+ const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2.0L);
/* tanh(x+iy) = (sinh(2x) + i*sin(2y))/(cosh(2x) + cos(2y))
= (sinh(x)*cosh(x) + i*sin(y)*cos(y))/(sinh(x)^2 + cos(y)^2). */
@@ -71,7 +73,7 @@ __ctanhl (__complex__ long double x)
the real part is +/- 1, the imaginary part is
sin(y)*cos(y)/sinh(x)^2 = 4*sin(y)*cos(y)/exp(2x). */
long double exp_2t = __ieee754_expl (2 * t);
- __real__ res = __copysignl (1.0, __real__ x);
+ __real__ res = __copysignl (1.0L, __real__ x);
__imag__ res = 4 * sinix * cosix;
__real__ x = fabsl (__real__ x);
__real__ x -= t;
@@ -83,22 +85,34 @@ __ctanhl (__complex__ long double x)
__imag__ res /= exp_2t;
}
else
- __imag__ res /= __ieee754_expl (2 * __real__ x);
+ __imag__ res /= __ieee754_expl (2.0L * __real__ x);
}
else
{
- long double sinhrx = __ieee754_sinhl (__real__ x);
- long double coshrx = __ieee754_coshl (__real__ x);
+ long double sinhrx, coshrx;
+ if (fabs (__real__ x) > LDBL_MIN)
+ {
+ sinhrx = __ieee754_sinhl (__real__ x);
+ coshrx = __ieee754_coshl (__real__ x);
+ }
+ else
+ {
+ sinhrx = __real__ x;
+ coshrx = 1.0L;
+ }
- den = sinhrx * sinhrx + cosix * cosix;
- __real__ res = sinhrx * coshrx / den;
- __imag__ res = sinix * cosix / den;
+ if (fabsl (sinhrx) > fabsl (cosix) * ldbl_eps)
+ den = sinhrx * sinhrx + cosix * cosix;
+ else
+ den = cosix * cosix;
+ __real__ res = sinhrx * (coshrx / den);
+ __imag__ res = sinix * (cosix / den);
}
/* __gcc_qmul does not respect -0.0 so we need the following fixup. */
- if ((__real__ res == 0.0) && (__real__ x == 0.0))
+ if ((__real__ res == 0.0L) && (__real__ x == 0.0L))
__real__ res = __real__ x;
- if ((__real__ res == 0.0) && (__imag__ x == 0.0))
+ if ((__real__ res == 0.0L) && (__imag__ x == 0.0L))
__imag__ res = __imag__ x;
}
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c b/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c
index 9d89bbe..1979f44 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c
@@ -25,6 +25,8 @@
#include <math_private.h>
+/* IBM long double GCC builtin sets LDBL_EPSILON == LDBL_DENORM_MIN */
+static const long double ldbl_eps = 1.232595164407831e-32; /* 2^-106 */
__complex__ long double
__ctanl (__complex__ long double x)
@@ -55,7 +57,7 @@ __ctanl (__complex__ long double x)
{
long double sinrx, cosrx;
long double den;
- const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2);
+ const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2.0L);
/* tan(x+iy) = (sin(2x) + i*sinh(2y))/(cos(2x) + cosh(2y))
= (sin(x)*cos(x) + i*sinh(y)*cosh(y)/(cos(x)^2 + sinh(y)^2). */
@@ -70,7 +72,7 @@ __ctanl (__complex__ long double x)
sin(x)*cos(x)/sinh(y)^2 = 4*sin(x)*cos(x)/exp(2y). */
long double exp_2t = __ieee754_expl (2 * t);
- __imag__ res = __copysignl (1.0, __imag__ x);
+ __imag__ res = __copysignl (1.0L, __imag__ x);
__real__ res = 4 * sinrx * cosrx;
__imag__ x = fabsl (__imag__ x);
__imag__ x -= t;
@@ -82,23 +84,35 @@ __ctanl (__complex__ long double x)
__real__ res /= exp_2t;
}
else
- __real__ res /= __ieee754_expl (2 * __imag__ x);
+ __real__ res /= __ieee754_expl (2.0L * __imag__ x);
}
else
{
- long double sinhix = __ieee754_sinhl (__imag__ x);
- long double coshix = __ieee754_coshl (__imag__ x);
+ long double sinhix, coshix;
+ if (fabsl (__imag__ x) > LDBL_MIN)
+ {
+ sinhix = __ieee754_sinhl (__imag__ x);
+ coshix = __ieee754_coshl (__imag__ x);
+ }
+ else
+ {
+ sinhix = __imag__ x;
+ coshix = 1.0L;
+ }
- den = cosrx * cosrx + sinhix * sinhix;
- __real__ res = sinrx * cosrx / den;
- __imag__ res = sinhix * coshix / den;
+ if (fabsl (sinhix) > fabsl (cosrx) * ldbl_eps)
+ den = cosrx * cosrx + sinhix * sinhix;
+ else
+ den = cosrx * cosrx;
+ __real__ res = sinrx * (cosrx / den);
+ __imag__ res = sinhix * (coshix / den);
}
/* __gcc_qmul does not respect -0.0 so we need the following fixup. */
- if ((__real__ res == 0.0) && (__real__ x == 0.0))
+ if ((__real__ res == 0.0L) && (__real__ x == 0.0L))
__real__ res = __real__ x;
- if ((__real__ res == 0.0) && (__imag__ x == 0.0))
+ if ((__real__ res == 0.0L) && (__imag__ x == 0.0L))
__imag__ res = __imag__ x;
}
diff --git a/sysdeps/powerpc/fpu/libm-test-ulps b/sysdeps/powerpc/fpu/libm-test-ulps
index 66576a5..ffb8e3a 100644
--- a/sysdeps/powerpc/fpu/libm-test-ulps
+++ b/sysdeps/powerpc/fpu/libm-test-ulps
@@ -1336,33 +1336,120 @@ ldouble: 1
Test "Real part of: ctan (0x1p1023 + 1 i) == -0.2254627924997545057926782581695274244229 + 0.8786063118883068695462540226219865087189 i":
double: 1
idouble: 1
+Test "Imaginary part of: ctan (0x1p1023 + 1 i) == -0.2254627924997545057926782581695274244229 + 0.8786063118883068695462540226219865087189 i":
+ildouble: 1
+ldouble: 1
Test "Real part of: ctan (0x1p127 + 1 i) == 0.2446359391192790896381501310437708987204 + 0.9101334047676183761532873794426475906201 i":
float: 1
ifloat: 1
+ildouble: 1
+ldouble: 1
Test "Imaginary part of: ctan (0x1p127 + 1 i) == 0.2446359391192790896381501310437708987204 + 0.9101334047676183761532873794426475906201 i":
double: 1
float: 1
idouble: 1
ifloat: 1
+ildouble: 2
+ldouble: 2
Test "Real part of: ctan (0x3.243f6cp-1 + 0 i) == -2.287733242885645987394874673945769518150e7 + 0.0 i":
float: 1
ifloat: 1
+ildouble: 2
+ldouble: 2
Test "Real part of: ctan (1 + 47 i) == 2.729321264492904590777293425576722354636e-41 + 1.0 i":
ildouble: 2
ldouble: 2
+# ctan_downward
+Test "Real part of: ctan_downward (0x1.921fb54442d18p+0 + 0x1p-1074 i) == 1.633123935319536975596773704152891653086e16 + 1.317719414943508315995636961402669067843e-291 i":
+ildouble: 3
+ldouble: 3
+Test "Real part of: ctan_downward (0x1.921fb6p+0 + 0x1p-149 i) == -2.287733242885645987394874673945769518150e7 + 7.334008549954377778731880988481078535821e-31 i":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 4
+ldouble: 4
+Test "Imaginary part of: ctan_downward (0x1.921fb6p+0 + 0x1p-149 i) == -2.287733242885645987394874673945769518150e7 + 7.334008549954377778731880988481078535821e-31 i":
+float: 1
+ifloat: 1
+ildouble: 10
+ldouble: 10
+
+# ctan_tonearest
+Test "Real part of: ctan_tonearest (0x1.921fb6p+0 + 0x1p-149 i) == -2.287733242885645987394874673945769518150e7 + 7.334008549954377778731880988481078535821e-31 i":
+float: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
+Test "Imaginary part of: ctan_tonearest (0x1.921fb6p+0 + 0x1p-149 i) == -2.287733242885645987394874673945769518150e7 + 7.334008549954377778731880988481078535821e-31 i":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+# ctan_towardzero
+Test "Real part of: ctan_towardzero (0x1.921fb54442d18p+0 + 0x1p-1074 i) == 1.633123935319536975596773704152891653086e16 + 1.317719414943508315995636961402669067843e-291 i":
+ldouble: 4
+ildouble: 4
+Test "Imaginary part of: ctan_towardzero (0x1.921fb54442d18p+0 + 0x1p-1074 i) == 1.633123935319536975596773704152891653086e16 + 1.317719414943508315995636961402669067843e-291 i":
+ildouble: 13
+ldouble: 13
+Test "Real part of: ctan_towardzero (0x1.921fb6p+0 + 0x1p-149 i) == -2.287733242885645987394874673945769518150e7 + 7.334008549954377778731880988481078535821e-31 i":
+float: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
+Test "Imaginary part of: ctan_towardzero (0x1.921fb6p+0 + 0x1p-149 i) == -2.287733242885645987394874673945769518150e7 + 7.334008549954377778731880988481078535821e-31 i":
+float: 1
+ifloat: 1
+ildouble: 10
+ldouble: 10
+
+# ctan_upward
+Test "Real part of: ctan_upward (0x1.921fb54442d18p+0 + 0x1p-1074 i) == 1.633123935319536975596773704152891653086e16 + 1.317719414943508315995636961402669067843e-291 i":
+double: 1
+idouble: 1
+ildouble: 6
+ldouble: 6
+Test "Imaginary part of: ctan_upward (0x1.921fb54442d18p+0 + 0x1p-1074 i) == 1.633123935319536975596773704152891653086e16 + 1.317719414943508315995636961402669067843e-291 i":
+ildouble: 10
+ldouble: 10
+Test "Real part of: ctan_upward (0x1.921fb6p+0 + 0x1p-149 i) == -2.287733242885645987394874673945769518150e7 + 7.334008549954377778731880988481078535821e-31 i":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 6
+ldouble: 6
+Test "Imaginary part of: ctan_upward (0x1.921fb6p+0 + 0x1p-149 i) == -2.287733242885645987394874673945769518150e7 + 7.334008549954377778731880988481078535821e-31 i":
+double: 1
+float: 2
+idouble: 1
+ifloat: 2
+ildouble: 1
+ldouble: 1
+
# ctanh
Test "Real part of: ctanh (-2 - 3 i) == -0.965385879022133124278480269394560686 + 0.988437503832249372031403430350121098e-2 i":
double: 1
float: 2
idouble: 1
ifloat: 2
+idouble: 2
+ildouble: 2
+ldouble: 2
Test "Imaginary part of: ctanh (-2 - 3 i) == -0.965385879022133124278480269394560686 + 0.988437503832249372031403430350121098e-2 i":
double: 1
idouble: 1
+ildouble: 2
+ldouble: 2
Test "Imaginary part of: ctanh (0 + 0x3.243f6cp-1 i) == 0.0 - 2.287733242885645987394874673945769518150e7 i":
float: 1
ifloat: 1
+ildouble: 2
+ldouble: 2
Test "Imaginary part of: ctanh (0 + pi/4 i) == 0.0 + 1.0 i":
double: 1
float: 1
@@ -1378,22 +1465,100 @@ float: 1
ifloat: 1
ildouble: 2
ldouble: 2
+Test "Real part of: ctanh (1 + 0x1p1023 i) == 0.8786063118883068695462540226219865087189 - 0.2254627924997545057926782581695274244229 i":
+ildouble: 1
+ldouble: 1
Test "Imaginary part of: ctanh (1 + 0x1p1023 i) == 0.8786063118883068695462540226219865087189 - 0.2254627924997545057926782581695274244229 i":
-double: 1
idouble: 1
+double: 1
Test "Real part of: ctanh (1 + 0x1p127 i) == 0.9101334047676183761532873794426475906201 + 0.2446359391192790896381501310437708987204 i":
double: 1
float: 1
idouble: 1
ifloat: 1
+ildouble: 2
+ldouble: 2
Test "Imaginary part of: ctanh (1 + 0x1p127 i) == 0.9101334047676183761532873794426475906201 + 0.2446359391192790896381501310437708987204 i":
double: 1
float: 1
ifloat: 1
+ildouble: 1
+ldouble: 1
Test "Imaginary part of: ctanh (47 + 1 i) == 1.0 + 2.729321264492904590777293425576722354636e-41 i":
ildouble: 2
ldouble: 2
+# ctanh_downward
+Test "Imaginary part of: ctanh_downward (0x1p-1074 + 0x1.921fb54442d18p+0 i) == 1.317719414943508315995636961402669067843e-291 + 1.633123935319536975596773704152891653086e16 i":
+ildouble: 3
+ldouble: 3
+Test "Real part of: ctanh_downward (0x1p-149 + 0x1.921fb6p+0 i) == 7.334008549954377778731880988481078535821e-31 - 2.287733242885645987394874673945769518150e7 i":
+float: 1
+ifloat: 1
+ildouble: 10
+ldouble: 10
+Test "Imaginary part of: ctanh_downward (0x1p-149 + 0x1.921fb6p+0 i) == 7.334008549954377778731880988481078535821e-31 - 2.287733242885645987394874673945769518150e7 i":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 4
+ldouble: 4
+
+# ctanh_tonearest
+Test "Real part of: ctanh_tonearest (0x1p-149 + 0x1.921fb6p+0 i) == 7.334008549954377778731880988481078535821e-31 - 2.287733242885645987394874673945769518150e7 i":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: ctanh_tonearest (0x1p-149 + 0x1.921fb6p+0 i) == 7.334008549954377778731880988481078535821e-31 - 2.287733242885645987394874673945769518150e7 i":
+float: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
+
+# ctanh_towardzero
+Test "Real part of: ctanh_towardzero (0x1p-1074 + 0x1.921fb54442d18p+0 i) == 1.317719414943508315995636961402669067843e-291 + 1.633123935319536975596773704152891653086e16 i":
+ildouble: 13
+ldouble: 13
+Test "Imaginary part of: ctanh_towardzero (0x1p-1074 + 0x1.921fb54442d18p+0 i) == 1.317719414943508315995636961402669067843e-291 + 1.633123935319536975596773704152891653086e16 i":
+ildouble: 4
+ldouble: 4
+Test "Real part of: ctanh_towardzero (0x1p-149 + 0x1.921fb6p+0 i) == 7.334008549954377778731880988481078535821e-31 - 2.287733242885645987394874673945769518150e7 i":
+float: 1
+ifloat: 1
+ildouble: 10
+ldouble: 10
+Test "Imaginary part of: ctanh_towardzero (0x1p-149 + 0x1.921fb6p+0 i) == 7.334008549954377778731880988481078535821e-31 - 2.287733242885645987394874673945769518150e7 i":
+float: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
+
+# ctanh_upward
+Test "Real part of: ctanh_upward (0x1p-1074 + 0x1.921fb54442d18p+0 i) == 1.317719414943508315995636961402669067843e-291 + 1.633123935319536975596773704152891653086e16 i":
+ildouble: 10
+ldouble: 10
+Test "Imaginary part of: ctanh_upward (0x1p-1074 + 0x1.921fb54442d18p+0 i) == 1.317719414943508315995636961402669067843e-291 + 1.633123935319536975596773704152891653086e16 i":
+double: 1
+idouble: 1
+ildouble: 6
+ldouble: 6
+Test "Real part of: ctanh_upward (0x1p-149 + 0x1.921fb6p+0 i) == 7.334008549954377778731880988481078535821e-31 - 2.287733242885645987394874673945769518150e7 i":
+double: 1
+float: 2
+idouble: 1
+ifloat: 2
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: ctanh_upward (0x1p-149 + 0x1.921fb6p+0 i) == 7.334008549954377778731880988481078535821e-31 - 2.287733242885645987394874673945769518150e7 i":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 3
+ldouble: 3
+
# erf
Test "erf (1.25) == 0.922900128256458230136523481197281140":
double: 1
@@ -2714,9 +2879,63 @@ double: 1
float: 1
idouble: 1
ifloat: 1
+ildouble: 2
+ldouble: 2
+
+Function: Real part of "ctan_downward":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 4
+ldouble: 4
+
+Function: Imaginary part of "ctan_downward":
+float: 1
+ifloat: 1
+ildouble: 10
+ldouble: 10
+
+Function: Real part of "ctan_tonearest":
+float: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
+
+Function: Imaginary part of "ctan_tonearest":
+float: 1
+ifloat: 1
ildouble: 1
ldouble: 1
+Function: Real part of "ctan_towardzero":
+float: 1
+ifloat: 1
+ildouble: 4
+ldouble: 4
+
+Function: Imaginary part of "ctan_towardzero":
+float: 1
+ifloat: 1
+ildouble: 13
+ldouble: 13
+
+Function: Real part of "ctan_upward":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 6
+ldouble: 6
+
+Function: Imaginary part of "ctan_upward":
+double: 1
+float: 2
+idouble: 1
+ifloat: 2
+ildouble: 10
+ldouble: 10
+
Function: Real part of "ctanh":
double: 1
float: 2
@@ -2733,6 +2952,60 @@ ifloat: 1
ildouble: 2
ldouble: 2
+Function: Real part of "ctanh_downward":
+float: 1
+ifloat: 1
+ildouble: 10
+ldouble: 10
+
+Function: Imaginary part of "ctanh_downward":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 4
+ldouble: 4
+
+Function: Real part of "ctanh_tonearest":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: Imaginary part of "ctanh_tonearest":
+float: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
+
+Function: Real part of "ctanh_towardzero":
+float: 1
+ifloat: 1
+ildouble: 13
+ldouble: 13
+
+Function: Imaginary part of "ctanh_towardzero":
+float: 1
+ifloat: 1
+ildouble: 4
+ldouble: 4
+
+Function: Real part of "ctanh_upward":
+double: 1
+float: 2
+idouble: 1
+ifloat: 2
+ildouble: 10
+ldouble: 10
+
+Function: Imaginary part of "ctanh_upward":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 6
+ldouble: 6
+
Function: "erf":
double: 1
idouble: 1
--
1.7.1
--
Adhemerval Zanella Netto
Software Engineer
Linux Technology Center Brazil
Toolchain / GLIBC on Power Architecture
azanella@linux.vnet.ibm.com / azanella@br.ibm.com
+55 61 8642-9890