lgamma_r and lgammaf_r return the wrong sign for -0.f
brde at optusnet.com.au
Sat Sep 13 06:00:33 UTC 2014
On Fri, 12 Sep 2014, enh via freebsd-numerics wrote:
> On Fri, Sep 12, 2014 at 3:05 PM, Steve Kargl
> <sgk at troutmask.apl.washington.edu> wrote:
>> On Fri, Sep 12, 2014 at 02:15:37PM -0700, enh via freebsd-numerics wrote:
>>> if I pass -0.f to lgammaf_r, the sign returned in *signgamp is 1. this
>>> is incorrect --- it should be -1.
>>> both lgamma_r and lgammaf_r are affected, but the other special cases
>>> in those functions look fine to me.
>>> this is fixed in OpenBSD and glibc, but FreeBSD and NetBSD both have
>>> the same bug.
>> Are you sure FreeBSD has a bug? From n1256.pdf (committee draft of C99)
C99 doesn't even have signgam. But it has the usual complications for
+-0 when the sign of the result can be determined for a real function,
so it has tgamma(+-0) = +-Inf. POSIX specifies that signgam has the
same sign as gamma(), where gamma is spelled with a greek letter. (Bug;
it should say the same sign as tgamma(). The one spelled with the greek
letter should mean the mathematical one, and that has a complex pole at
0. Complex poles don't have signs.)
I looked at what tgamma() does. It is correct for +-0, but for large
finite args and +Inf it returns 1/zero. This is only correct for +Inf.
>>> patch below (whitespace mangled courtesy of gmail). i'd prefer to wait
>>> for this to be fixed in FreeBSD and pull down the fix rather than just
>>> fix it locally.
The patch can possibly be improved a little by duplicating the zero check
((ix|lx) == 0) instead of duplicating the sign check (hx < 0). This
gives the same number of checks, but 1 fewer setting of signgam.
In logl(), the checks are arranged more subtly to try to optimize for
the usual case. This works best in ld80 long double precision. The
single check (hx <= 0) classifies all negative and 0 x. This saves
as much as 1 whole cycle or 1.5% on modern x86. It works best in for
ld80 since hx doesn't contain any mantissa bits then. When hx contains
mantissa bits, its sign bit is usually checked using the C '<' operator
as in (hx < 0) and not as a bit. Sometimes the compiler can combine
this with another operator and do it as a bit test. The ((ix|lx) ==
0) check may be too special for this to be possible except in 64-bit
mode where hx, ix, and lx can be treated as a single 64-bit value (then
(hlx = 0) works, where int64_t hlx is hx and lx combined). The
optimization is going the other way an converting several bit tests
into a single signed comparison.
More information about the freebsd-numerics