This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix e_logl (128ibm) spurious underflow
- From: Adhemerval Zanella <azanella at linux dot vnet dot ibm dot com>
- To: "GNU C. Library" <libc-alpha at sourceware dot org>
- Date: Wed, 27 Feb 2013 17:11:23 -0300
- Subject: Fix e_logl (128ibm) spurious underflow
This patch fixes the spurious underflow exceptions triggered by cacos,
casin, and casinh test for IBM long double. This is in fact an issue with
__ieee754_logl implementation: in the domain reduction phase, where
the algorithm transform log(u) -> log(t) + log(u/t) and log(u/t) -> log(1+(u-t)/t),
if the input is high enough the calculation '1+(u-t)/t' might generate
subnormal numbers.
The patch I propose tests if the resulting coefficient used in the series
expansion might trigger a subnormal multiplication and if it is the case,
just set the expansion to 0. Since it is denormal multiplications and the
expect result is orders of magnitude higher than the series expansion, I
believe the loss of precision is acceptable.
I also added the ulps update for remaining cacos tests.
Any tips, comments, advices?
--
2013-02-27 Adhemerval Zanella <azanella@linux.vnet.ibm.com>
* sysdeps/ieee754/ldbl-128ibm/e_logl.c (__ieee754_logl): Fix spurious
underflow.
* sysdeps/powerpc/fpu/libm-test-ulps: Update.
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_logl.c b/sysdeps/ieee754/ldbl-128ibm/e_logl.c
index 14f47eb..15b5edf 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_logl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_logl.c
@@ -182,6 +182,9 @@ static const long double
ln2a = 6.93145751953125e-1L,
ln2b = 1.4286068203094172321214581765680755001344E-6L;
+static const long double
+ ldbl_epsilon = 0x1p-106L;
+
long double
__ieee754_logl(long double x)
{
@@ -258,7 +261,12 @@ __ieee754_logl(long double x)
}
/* Series expansion of log(1+z). */
w = z * z;
- y = ((((((((((((l15 * z
+ /* Avoid spurious underflows. */
+ if (__glibc_unlikely(w <= ldbl_epsilon))
+ y = 0.0L;
+ else
+ {
+ y = ((((((((((((l15 * z
+ l14) * z
+ l13) * z
+ l12) * z
@@ -271,7 +279,8 @@ __ieee754_logl(long double x)
+ l5) * z
+ l4) * z
+ l3) * z * w;
- y -= 0.5 * w;
+ y -= 0.5 * w;
+ }
y += e * ln2b; /* Base 2 exponent offset times ln(2). */
y += z;
y += logtbl[k-26]; /* log(t) - (t-1) */
diff --git a/sysdeps/powerpc/fpu/libm-test-ulps b/sysdeps/powerpc/fpu/libm-test-ulps
index 4221967..5072190 100644
--- a/sysdeps/powerpc/fpu/libm-test-ulps
+++ b/sysdeps/powerpc/fpu/libm-test-ulps
@@ -243,9 +243,11 @@ idouble: 1
Test "Real part of: cacos (-0.5 + +0 i) == 2.094395102393195492308428922186335256131 - 0 i":
double: 1
idouble: 1
+ldouble: 1
Test "Real part of: cacos (-0.5 - 0 i) == 2.094395102393195492308428922186335256131 + +0 i":
double: 1
idouble: 1
+ldouble: 1
Test "Imaginary part of: cacos (-1.5 + +0 i) == pi - 0.9624236501192068949955178268487368462704 i":
double: 1
float: 1
@@ -265,6 +267,10 @@ double: 1
float: 1
idouble: 1
ifloat: 1
+Test "Imaginary part of: cacos (0x1.fp1023 + 0x1.fp1023 i) == 7.853981633974483096156608458198757210493e-1 - 7.107906849659093345062145442726115449315e2 i":
+ldouble: 1
+Test "Imaginary part of: cacos (0x1.fp127 + 0x1.fp127 i) == 7.853981633974483096156608458198757210493e-1 - 8.973081118419833726837456344608533993585e1 i":
+ldouble: 1
# cacosh
Test "Real part of: cacosh (+0 + 0.5 i) == 0.4812118250596034474977589134243684231352 + pi/2 i":
--
1.7.1