Date: Sat, 13 Sep 2014 16:00:22 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: enh <enh@google.com> Cc: freebsd-numerics@freebsd.org, Steve Kargl <sgk@troutmask.apl.washington.edu> Subject: Re: lgamma_r and lgammaf_r return the wrong sign for -0.f Message-ID: <20140913152505.E950@besplex.bde.org> In-Reply-To: <CAJgzZoruuVxLJ1sdue-NmvxKJwW7WKrY6GsopKMQi3Us4goxCQ@mail.gmail.com> References: <CAJgzZopa-d=eR7zkqhffsjMY0NEavhqDA-B3V9bQdaJd6BMO2A@mail.gmail.com> <20140912220544.GA96714@troutmask.apl.washington.edu> <CAJgzZoruuVxLJ1sdue-NmvxKJwW7WKrY6GsopKMQi3Us4goxCQ@mail.gmail.com>
next in thread | previous in thread | raw e-mail | index | archive | help
On Fri, 12 Sep 2014, enh via freebsd-numerics wrote: > On Fri, Sep 12, 2014 at 3:05 PM, Steve Kargl > <sgk@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. Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20140913152505.E950>