This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: [PATCH] Update sparc ULPs.
- From: David Miller <davem at davemloft dot net>
- To: joseph at codesourcery dot com
- Cc: libc-alpha at sourceware dot org
- Date: Fri, 16 Nov 2012 01:21:26 -0500 (EST)
- Subject: Re: [PATCH] Update sparc ULPs.
- References: <Pine.LNX.4.64.1209241254140.23802@digraph.polyomino.org.uk><20120924.135523.2105566314867619791.davem@davemloft.net><Pine.LNX.4.64.1209241940230.31396@digraph.polyomino.org.uk>
From: "Joseph S. Myers" <joseph@codesourcery.com>
Date: Mon, 24 Sep 2012 19:48:44 +0000
> If x is small enough (absolute value less than LDBL_EPSILON / 2.0L, I
> think) then the correct rounded-to-nearest result of log1pl (x) is just x,
> so code can test for that case and avoid the more complicated calculations
> with potential underflow (while taking care to generate underflow if
> correct). Properly, that's x with inexact exception (if not 0), and
> underflow exception if x is subnormal (or, on tinyness-before-rounding
> machines, if x is +LDBL_MIN, but that detail isn't something I'd expect
> glibc to implement).
Indeed, that seems to do the trick. I used an int conversion to
force the exception as that seems to be the cheapest way to do
this in the soft-fp library.
What I really don't understand is why ldbl-128ibm's variant of this
code doesn't trigger this same problem on powerpc testsuite runs.
Either it is happening (and therefore ldbl-128ibm/s_log1pl.c needs a
similar change to this) or powerpc isn't generating the exceptions.
Anyways, any objections to this patch?
2012-11-15 David S. Miller <davem@davemloft.net>
* sysdeps/ieee754/ldbl-128/s_log1pl.c (__log1pl): If xm1 is
smaller than LDBL_EPSILON/2.0L, just return xm1.
diff --git a/sysdeps/ieee754/ldbl-128/s_log1pl.c b/sysdeps/ieee754/ldbl-128/s_log1pl.c
index 4ecea0f..fad18e9 100644
--- a/sysdeps/ieee754/ldbl-128/s_log1pl.c
+++ b/sysdeps/ieee754/ldbl-128/s_log1pl.c
@@ -138,6 +138,12 @@ __log1pl (long double xm1)
&& (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
return xm1;
+ if ((hx & 0x7fffffff) < 0x3f8e0000)
+ {
+ if ((int) xm1 == 0)
+ return xm1;
+ }
+
x = xm1 + 1.0L;
/* log1p(-1) = -inf */