diff mbox

[BZ,#16823] Fix log1pl returning wrong infinity sign

Message ID 53562754.3040706@linux.vnet.ibm.com
State Committed
Headers

Commit Message

Stefan Liebler April 22, 2014, 8:24 a.m. UTC 
  Hi,

log1pl (-1) now returns with -inf instead of +inf in FE_DOWNWARD 
rounding mode. The source code producing the -inf for the double/float 
variants is the same, but on S390 the compiler generates a division by 
constant +zero and ignores the (x-x).
Thus i also changed it to the constant zero in the double/float functions.
I´ve tested it on S390.
Please test it on other platforms because there are also reports for the 
same failure on aarch64 (see bug comment 1)
and sparc (https://www.sourceware.org/ml/libc-alpha/2014-04/msg00410.html).

Give feedback in order to commit the patch.

Bye

---
2014-04-22  Stefan Liebler  <stli@linux.vnet.ibm.com>

	[BZ #16823]
	* sysdeps/ieee754/ldbl-128/s_log1pl.c (__log1pl):
	Always divide by positive zero when computing -Inf result.
	* sysdeps/ieee754/dbl-64/s_log1p.c (__log1p): Likewise.
	* sysdeps/ieee754/flt-32/s_log1pf.c (__log1pf): Likewise.
---

On 04/09/2014 03:24 PM, Stefan Liebler wrote:
> Hi,
>
> on S390 test-double fails for log1pl(-1) with infinity has wrong sign in
> rounding mode FE_DOWNWARD. See Bug 16823.
> In this rounding mode, (x-x) = -0.
> In all other rounding modes, (x-x) = +0.
> A division with divisor -0 leads to +inf, while +0 results in -inf as
> expected. This patch changes the divisor to a const +0.
> Tested on s390/s390x.
>
> Is this okay?
>
> Bye
>
> ---
> 2014-04-09  Stefan Liebler  <stli@linux.vnet.ibm.com>
>
>      [BZ #16823]
>      * sysdeps/ieee754/ldbl-128/s_log1pl.c (__log1pl):
>      Use const positive zero instead of (x-x).
> ---

Comments

David Miller April 23, 2014, 12:42 a.m. UTC | #1 
From: Stefan Liebler <stli@linux.vnet.ibm.com>
Date: Tue, 22 Apr 2014 10:24:52 +0200

> log1pl (-1) now returns with -inf instead of +inf in FE_DOWNWARD
> rounding mode. The source code producing the -inf for the double/float
> variants is the same, but on S390 the compiler generates a division by
> constant +zero and ignores the (x-x).
> Thus i also changed it to the constant zero in the double/float
> functions.
> I´ve tested it on S390.
> Please test it on other platforms because there are also reports for
> the same failure on aarch64 (see bug comment 1)
> and sparc
> (https://www.sourceware.org/ml/libc-alpha/2014-04/msg00410.html).
> 
> Give feedback in order to commit the patch.

I can confirm that this fixes the testsuite failures on sparc.
Marcus Shawcroft April 28, 2014, 9:02 a.m. UTC | #2 
On 22 April 2014 09:24, Stefan Liebler <stli@linux.vnet.ibm.com> wrote:
> Hi,
>
> log1pl (-1) now returns with -inf instead of +inf in FE_DOWNWARD rounding
> mode. The source code producing the -inf for the double/float variants is
> the same, but on S390 the compiler generates a division by constant +zero
> and ignores the (x-x).
> Thus i also changed it to the constant zero in the double/float functions.
> I´ve tested it on S390.
> Please test it on other platforms because there are also reports for the
> same failure on aarch64 (see bug comment 1)
> and sparc (https://www.sourceware.org/ml/libc-alpha/2014-04/msg00410.html).
>
> Give feedback in order to commit the patch.

This resolves the failure on aarch64....
Cheers
/Marcus
Andreas Krebbel April 29, 2014, 1:46 p.m. UTC | #3 
On 04/22/2014 10:24 AM, Stefan Liebler wrote:
> 2014-04-22  Stefan Liebler  <stli@linux.vnet.ibm.com>
> 
> 	[BZ #16823]
> 	* sysdeps/ieee754/ldbl-128/s_log1pl.c (__log1pl):
> 	Always divide by positive zero when computing -Inf result.
> 	* sysdeps/ieee754/dbl-64/s_log1p.c (__log1p): Likewise.
> 	* sysdeps/ieee754/flt-32/s_log1pf.c (__log1pf): Likewise.

Applied.  Thanks!

-Andreas-
diff mbox

Patch

diff --git a/sysdeps/ieee754/dbl-64/s_log1p.c b/sysdeps/ieee754/dbl-64/s_log1p.c
index fd4dce5..c922148 100644
--- a/sysdeps/ieee754/dbl-64/s_log1p.c
+++ b/sysdeps/ieee754/dbl-64/s_log1p.c
@@ -110,7 +110,7 @@  __log1p (double x)
       if (__glibc_unlikely (ax >= 0x3ff00000))           /* x <= -1.0 */
 	{
 	  if (x == -1.0)
-	    return -two54 / (x - x);            /* log1p(-1)=+inf */
+	    return -two54 / zero;               /* log1p(-1)=-inf */
 	  else
 	    return (x - x) / (x - x);           /* log1p(x<-1)=NaN */
 	}
diff --git a/sysdeps/ieee754/flt-32/s_log1pf.c b/sysdeps/ieee754/flt-32/s_log1pf.c
index 0307277..5f00feb 100644
--- a/sysdeps/ieee754/flt-32/s_log1pf.c
+++ b/sysdeps/ieee754/flt-32/s_log1pf.c
@@ -42,7 +42,7 @@  __log1pf(float x)
 	k = 1;
 	if (hx < 0x3ed413d7) {			/* x < 0.41422  */
 	    if(ax>=0x3f800000) {		/* x <= -1.0 */
-		if(x==(float)-1.0) return -two25/(x-x); /* log1p(-1)=+inf */
+		if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=-inf */
 		else return (x-x)/(x-x);	/* log1p(x<-1)=NaN */
 	    }
 	    if(ax<0x31000000) {			/* |x| < 2**-29 */
diff --git a/sysdeps/ieee754/ldbl-128/s_log1pl.c b/sysdeps/ieee754/ldbl-128/s_log1pl.c
index d991e8a..d8d89f0 100644
--- a/sysdeps/ieee754/ldbl-128/s_log1pl.c
+++ b/sysdeps/ieee754/ldbl-128/s_log1pl.c
@@ -150,7 +150,7 @@  __log1pl (long double xm1)
   if (x <= 0.0L)
     {
       if (x == 0.0L)
-	return (-1.0L / (x - x));
+	return (-1.0L / zero);  /* log1p(-1) = -inf */
       else
 	return (zero / (x - x));
     }